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(%);

@leafgreen This is just a comment. Try the following:
 

restart; 
u:=-(x[1]-x[2])/((x[1]-x[2])^2+(y[1]-y[2])^2)-(x[1]-x[3])/((x[1]-x[3])^2+(y[1]-y[3])^2);
F:=unapply(u,x[1], y[1], x[2], y[2], x[3], y[3]);
f:=(x,y)->-(x[1]-x[2])/((x[1]-x[2])^2+(y[1]-y[2])^2)-(x[1]-x[3])/((x[1]-x[3])^2+(y[1]-y[3])^2);

You will see that unapply uses other symbols in defining the function. But you get the function you want.
In the second example (f) the function f is a function of x and y (typically vectors). This syntax is used e.g. by dsolve/numeric behind the scene, but can also be used directly by the user.

You should upload the worksheet int the MaplePrimes editor. Use the fat, green arrow to do that.

@ernilesh80 You have equality signs ( = ) instead of assignments (:= ) in the line that should assign to best_value.

Incidentally, you have defined best_value as a list, not as an array. You are allowed to assign to small lists (length less than 101), but it is a bad idea since an entirely new list is created every time. Lists are not mutable as opposed to arrays, matrices, vectors, and tables.
So use an array instead (here basically for matters of principle).
best_value:=Array([0,0,0,0,0]);
### Try this:
 

restart;
L:=[0,0,0,0,0];
L[1]:=1;
L;
L:=[0$101];
L[1]:=1;

@ernilesh80 You forgot to do
with(LinearAlgebra);
before the loop, e.g. at the very top.
On the other hand I would use the long form for the Hessian VectorCalculus:-Hessian instead of using with(VectorCalculus). That package redefines all kinds of normal things, which I find disturbing.

@Harry Garst This clearly is a bug in minimize.
I believe that it can be traced from minimize to
`minimize/cell/check`, from there to
`minimize/cell/internal`, then to
`minimize/cell/decouple` then to
`minimize/cell/sum` then to
`minimize/cell/variate` then to
`minimize/cell/univariate`
after which I gave up.
###################################################
A few observations:
1. If you copy the final temp2 and paste it in a fresh worksheet then there is no problem:
 

temp2:=(-3.79741609404203+1.21140500697585*sqrt(v^2))^2+(-1.33163418865668+1.21140500697585*sqrt((1-v)^2))^2+(-.865852283271325+1.21140500697585*sqrt((2-v)^2))^2+(-1.40007037788597+1.21140500697585*sqrt((3-v)^2))^2+(-.934288472500622+1.21140500697585*sqrt((4-v)^2))^2+(-2.46850656711527+1.21140500697585*sqrt((5-v)^2))^2+(-4.00272466172992+1.21140500697585*sqrt((6-v)^2))^2+(-4.53694275634457+1.21140500697585*sqrt((7-v)^2))^2+(-7.07116085095921+1.21140500697585*sqrt((8-v)^2))^2+(-7.60537894557386+1.21140500697585*sqrt((9-v)^2))^2;
minimize(temp2,v=1..8,location);

2. If you use abs(xxx) instead of sqrt(xxx^2) in the model as well as in temp2 then there is no problem.
3. During my wild goose chase through subprocedures of minimize I saw the name forget_Proc come up. This made me try this:
 

restart;
X := Vector([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]):
d := 2.730469991:
Y:=Vector( [5, 9, 11, 12, 14, 14, 14, 15, 14, 15]):
m := Statistics:-LinearFit(a+b*t+c*sqrt((t-d)^2), X, Y, t);
#m2 := Statistics:-LinearFit(a+b*t+c*abs(t-d), X, Y, t);
forget(minimize);

The response from forget minimize was the error message:

Error, (in forget) lexically scoped parameter out of context
If you uncomment the m2 line, comment the m line, and start over, then there is no error from forget.

I shall submit an SCR.

Could you provide us with all the code in the form of text or an uploaded worksheet?
Images don't help much.
Nobody is going to type the stuff and they won't even know how temp2 was defined or what happened before.