@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.

@ernilesh80 Since you don't know how large the two arrays have to be you can just define them as arrays:
pd_arr:=Array();
id_arr:=Array();
###
Then in the assignments to pd_arr and id_arr in the loops you must use parentheses instead of square brackets as in
pd_arr(pd, 1) := op_value[1]; pd_arr(pd, 2) := op_value[2]; pd_arr(pd, 3) := op_value[3]; pd_arr(pd, 4) := op_value[4];
Do similarly for id_arr.
The loops take a while since
nops~([T_arr,E_arr,W_arr,p_arr]);
`*`(op(%)); 
returns 7436529.
I stopped the computation after a while and inspected both arrays by just doing:
pd_array;  # at that time a 1..8667x1..4 array
id_array ;  # at that time a 1..46980x1..4 array, i.e. larger than the 1..10000x1..4 you used.
## Incidentally, it may not be more efficient to do this, but it is shorter:
pd_arr(pd, 1..4) := Array(op_value);
and similarly for id_arr.

@memdream You should be aware that (-1)^(1/4) in Maple is the principal root:
evalc( (-1)^(1/4) );
                                 (1/2)*sqrt(2)+(1/2*I)*sqrt(2)

@Rouben Rostamian  Yes, when I paste the 2D input on a 1D input line I get
exp(4.605170186*`9`)
The factor `9` is understood by Maple as a name, not as the integer 9.