@Christopher2222 You could use the following script in your maple.ini file:
 

ModGasDensity:=proc(elem) local vuu;
  uses ScientificConstants;
  if not HasElement(elem) then error "Element %1 is not included", elem end if;
  vuu:=rhs(GetElement(elem,'density')[2]);
  if convert(subs(vuu,units),`global`) = 'kg/L' then
    vuu:=[vuu[1],vuu[2],'units'='g/L']; 
    printf("Modifying the unit for the density of the gas %a from kg/L to g/L\n",elem);
    ModifyElement(elem,'density'=vuu)
  else
    NULL
  end if
end proc:
for i from 1 to 110 do 
   try
     ModGasDensity(i)
   catch:
   end try
end do;

You may want to comment out the printing.
Test:

with(ScientificConstants);
GetValue(Element(Krypton,density(gas)));
Units:-UsingSystem();

### Since Maple 8 through 15 have L where Maple 16 and up have kg/L I modified the script some and included it in a procedure call to avoid assignments to global variables.
This should work in all versions (I hope). At least it works in Maple 12, 15, 16, 2017.
In Maple 8 we cannot use 'uses ScientificConstants' so the long form have to be used for HasElement and the rest. I have done that, but another problem with 8 is that it doesn't read the maple.ini file in my Users folder. Not a big problem since I only look at Maple 8 for purposes of comparison with later versions.

proc() local VERSION,ns,N, ModGasDensity,i;
   VERSION:=convert(kernelopts(version),string):
   ns:=StringTools:-Search(",",VERSION):
   N:=parse(VERSION[7..ns-1]);
   if N<8 then return NULL end if;
   ModGasDensity:=proc(elem) local vuu,U;
     if not ScientificConstants:-HasElement(elem) then error "Element %1 is not included", elem end if;
     vuu:=rhs(ScientificConstants:-GetElement(elem,'density')[2]);
     U:=convert(subs(vuu,units),`global`);
     if N>=16 and U = 'kg/L' or N<16 and U='L' then
       vuu:=[vuu[1],vuu[2],'units'='g/L'];
       printf("Modifying the unit for the density of the gas %a to g/L\n",elem);
       ScientificConstants:-ModifyElement(elem,'density'=vuu)
     else
       NULL
     end if
   end proc;
   for i from 1 to 120 do
     try
       ModGasDensity(i)
     catch:
     end try
   end do;
end proc();

@Mac Dude Right, it ought to be easy to fix.
To see the units used in the record it is easier to do
ScientificConstants:-GetElement(H,density);
    

In Maple 8, 12 and 15 (probably 9, 9.5, 10, 11, 13, and 14 too) the last item is units = L.
So yes, certainly the part of ScientificConstants having to do with the density of gases couldn't have been used much over the years: surviving from Maple 8 (when the package was introduced) to Maple 15 with a density reported as 1000 times too low, and then from Maple 16 to Maple 2017 with a density that is 1000 times too large (both when SI units are desired).
 

A comment about the unit used at the time the data were put into the system:
In Maple 15 the unit reported by line 42 in the procedure GetValue  by
ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units
is L, i.e. liter, a measure of volume! The 0.088L is then correctly converted to the SI unit m^3 as 0.000088.
##
In Maple 16 (and following versions) this should have been corrected to the unit g/L, but was made kg/L, resulting in 88.0.

## Note x g/L = x kg/m^3.

@9009134 You have a problem with a wide range in the numerical coefficients:
Take a look at
 

(I set Digits:=15)
sys3;
SYS;
## The absolute values of the floats in the systems:
abs~(indets(evalf(sys3),float));
(min,max)(%); # Answer: 4.50605975510204*10^(-80), 4.66666666666667*10^7
abs~(indets(evalf(SYS),float));
(min,max)(%); # Answer 1.32275132275132*10^(-41), 4.82434787943840*10^99

I don't know if you will be able to scale this problem so that this huge relative difference in size disappears or at least gets considerably smaller. Otherwise you will have to work with a large value of Digits, way into the software floating point range, i.e. much larger than 15, which will slow down the handling considerably.

@9009134 You must make sure that the different quantities actually get defined before embarking on the dsolve process.
When I execute your worksheet ss.mw by using !!! in the top panel I find that SYS evaluates to NULL, i.e. the solve command has no solution. Furthermore newsys2 evaluates to

which is not a system of odes.
So my suggestion is to write the whole thing in a worksheet (not a document) using 1D math input (aka Maple input). 
Write one command at a time and examine the output to see if there is anything wrong at each step.
As it is right now it is a mess (quite frankly).

@Dmitrii It is highly unlikely that you would want a function just defined on your list x having corresponding values in your list y. But just in case you do, that is easily done:
 

restart;
x:=[0,1,2,3,4,5];
y:=[0,2,4,6,8,10];
for i from 1 to numelems(x) do f(x[i]):=y[i] end do;
eval(f); #f has option remember
op(4,eval(f)); # The remember table
f(t); # returns unevaluated.
f(0.5); # returns unevaluated since f(0.5) was never defined.
f(2); # returns 4 since 2 is in the remember table.

Surely, you want some kind of interpolation as Rouben is saying.

If you are pretty sure that this behavior applies to any Maple product I would very much like to see the example in Maple itself.
Give us text or upload a worksheet, please.

@John Fredsted Maybe Tom Leslie's example was too simple to show the problem clearly.
Anyway, I made a slightly more complicated example, where the first method of solving successively is fine, but the second is not and should be avoided:
 

restart;
ode1:=diff(x(t),t)=x(t);
res1:=dsolve(ode1);
ode2:=diff(y(t),t)=y(t)+x(t);
res2:=dsolve(eval(ode2,res1)); #Fine
## And is the same as solved as a system:
dsolve({ode1,ode2}); 
########################
## This way is no good:
dsolve(ode2,y(t)); 
eval(%,res1);
value(%);

@tomleslie I wouldn't complain about your first example either. It amounts to
 

restart;
de1:=diff(y(x),x)=1;
de2:=diff(u(x),x)=1
dsolve(de1);
dsolve(de2);

In both cases the arbitrary constant is called _C1. That is well known and can be counted on.
Consistency is important. In other parts of Maple constants are renumbered every time you redo the calculation as in _Z3 and next time _Z4 etc. which is rather annoying.

to 10 do solve(sin(x)=0,allsolutions) end do;

@Pascal4QM Since overloading `+` breaks existing behavior of `+` with lists of equal length:
restart;
[a, b, c]+[d, e, f];
                             [d+a, e+b, f+c]
it would be better to stick with Fabio's own suggestion and overload `union`:
 

restart;
ListOperations := module () export `union`; option package; 
  `union` := proc (a::list, b::list) option overload; [op(a), op(b)] end proc:
end module:
with(ListOperations);
[a, b, c] union [d, e, f];
[a, b, c]+[d, e, f]; #Still works as before overloading `union`

Your question is rather vague to me. Of course the mapping function in e.g.
conformal(z^2, z=0..2+2*I);
is just z->z^2.
So could you be a bit more specific?

Whereas I haven't been able to find anything in the code for invlaplace that makes any actual use of the fourth argument opt (if present), such is not the case for laplace.
inttrans[laplace] calls `inttrans/laplace/internal`, which in turn calls `inttrans/laplace/main`. The option opt is brought along and here is actual code handling the option NO_INT.
showstat(`inttrans/laplace/main`,29..35);

The corresponding train of events for invlaplace is:
From inttrans[invlaplace] to `inttrans/invlaplace/internal` to `inttrans/invlaplace/main`. Thus quite similar to laplace namewise.
But I cannot locate any code that actually handles the fourth argument opt although it is taken along for the ride.
Clearly I may have overlooked something.

@umar khan I'm puzzled by your epsilon(rho), which is defined piecewise.
With rho given as a fixed concrete value (as it is in your worksheet, rho=2.727373039e14 ), epsilon(rho) is just an equally fixed concrete value, in this case epsilon(rho) = 2.424644682*10^35.
##
So do you mean: How does the solution depend on the parameter E := epsilon(rho)?
If so you could replace epsilon(rho) in your odes by the name E and use the parameters option in dsolve/numeric.
You can combine the parameters option with the events option in the form used by Tom Leslie:

sol := dsolve({odes} union ics, numeric, parameters=[E],events=[[P(r),halt]]);
sol(parameters=[2.424644682*10^35]); #Testing the value you used before
plots:-odeplot(sol, [r, P(r)], r = 0.1e-5 .. 0.15e7);
## Making a procedure that sets the parameter value:
Q:=proc(E) sol(parameters=[E]); sol end proc;
## Animation in E:
plots:-animate(plots:-odeplot,['Q(E)',[r,P(r)],ep..1.11e7],E=10^34..3*10^34);

But probably you meant something entirely different?
 

Notice that each time you call foo a new local is created:
 

foo:=proc(n)
   local x;
   x^n;
end proc;
x1:=foo(1);
x2:=foo(1);
x1-x2;

x1-x2 prints as x-x and doesn't simplify to zero.
If you look at their addresses you will see that they are different:
addressof(x1);
addressof(x2);

@nm Clearly the problem on the interval -1..1 can be transformed to a problem on 0..2, which I shall do below.
I agree with Rouben that your argument in the pdf file fails when you assert that B must be taking to be zero since otherwise you have two unknowns A and B. But clearly you have infinitely many, since A (and B) must depend on lambda.
From your two equations (3) and (4) together as a linear system in the unknowns A and B follows (if you want a nontrivial solution) that the determinant is zero. The determinant is 2*sin(sqrt(lambda))*cos(sqrt(lambda)), which is sin(2*sqrt(lambda)).
Thus 2*sqrt(lambda) = n*Pi, where n is a positive integer. We still must use one of the two equations (3) or (4).
If you continue along this line you will end up with sine contributions as well as cosine contributions (in general).
## I give the code for the solution using transformation from -1..1 to 0..2.
 

restart;
## First we just document that solving on 0..L is done painlessly:
pde:=diff(u(x,t),t)=diff(u(x,t),x,x);
ibcs:=u(0,t)=0,u(L,t)=0,u(x,0)=f(x);
pdsolve({pde,ibcs}) assuming L>0;
## Now the present problem
icbs2:=u(-1,t)=0,u(1,t)=0,u(x,0)=f(x);
##Changing from -1..1 to 0..2:
sys:=PDEtools:-dchange({u(x,t)=v(xi,t),x=xi-1,f(x)=g(xi)},{pde,icbs2},[v,g,xi]);
res:=pdsolve(sys); # Done as in the beginning.
## Now back to the original variables:
evalindets(res,Integral,s->IntegrationTools:-Change(s,xi=y+1));
res:=u(x,t)=rhs(eval(%,{xi=x+1,g=(xi->f(xi-1))}));
N:=op(indets(res,suffixed(_Z,integer)));
## First taking the contribution from even N:
subs(N=2*m,rhs(res)); 
res1:=subsop(2=(m=1..infinity),%);
evalindets(res1,specfunc(sin),expand);
res1:=simplify(%) assuming m::posint; # redefining res1
## Now taking the contribution from odd N:
subs(N=2*m-1,rhs(res)); 
res2:=subsop(2=(m=1..infinity),%);
evalindets(res2,specfunc(sin),expand,2*m-1);
res2:=simplify(%) assuming m::posint; # redefining res2
## The final solution:
sol:=u(x,t)=res1+res2;
## Examples
eval(sol,f=(x->1-x^2)); #An even function
value(%);
eval(sol,f=(x->(1-x^2)*x)); #An odd function
value(%);
eval(sol,f=(x->(1-x^2)*(x+1))); #The sum of the two: neither even nor odd
value(%);