@das_goon OK, try this instead:

sys1 := [diff(c(x, t), t) = gDiffusion*10^5*diff(c(x, t), x$2), diff(zeta(x, t), t) = KDiffusion*10^6*diff(zeta(x, t), x$2)];
pds1:=pdsolve(eval(sys1[2],KDiffusion=1),{zeta(0, t) = .4, zeta(x, 0) = .4, D[1](zeta)(0, t) = 0},numeric,time=t,range=0..2000,spacestep=3);
res1:=pds1:-value(output=listprocedure);
Z:=subs(res1,zeta(x,t));
Z(2000,1);
Z(2000,t); #Returns unevaluated as it should
pds2:=pdsolve(eval(sys1[1],gDiffusion=1),{c(2000, t) = Z(2000,t), c(x, 0) = 0, D[1](c)(3000, t) = sin((1/100)*t)},numeric,time=t,spacestep=3);
pds1:-plot3d(t=0..1000); p1:=%:
pds2:-plot3d(t=0..1000); p2:=%:
plots:-display(p1,p2);

It seems that you are just asking if Maple have ways of finding the maximum of a function of a single variable symbolically. Maple can do the usual calculus needed to find a maximum of f(x). Find f'(x), solve f'(x)=0. Examine f''(x). etc. That there are parameters in the function is basically irrelevant. 

@samiyare Yes, I must have made a mistake.
Try this for a start from NBT=2 in the loop in the previous worksheet:
GUESS:=[T(eta) = -eta+.7*eta^2, u(eta) = 0.1e-2, phi(eta) = 0.09, phi0 = 0.09, p2 = 0.4, U(eta) = 0.001*eta, UF(eta) = 0.0001*eta];
LOOP:
res:=GUESS:
for NBT2 in [seq(2-.1*i,i=0..19)] do
  res := dsolve(subs(NBT=NBT2,{ueq2,Teq,eq3,eq4,eq5,U(0)=0,UF(0)=0,((-cbf*rhobf+cp*rhop)*UF(1)+ rhobf*cbf*U(1))/10000=p2,UF(1)=Phiavg*U(1),u(0)=lambda*D(u)(0),u(1)=-lambda*D(u)(1),D(T)(0)=-1,phi(0)=phi0,T(0)=0}), numeric,method=bvp[midrich],approxsoln=res);
  RES[NBT2]:=eval(res)
end do:

I couldn't download the file you attached, but used my own with the new parameters.
You may ask how I got GUESS this time.  I had success with my own unsophisticated code with some wild guesses and constructed GUESS based on that. My code is unsophisticated in the sense that it reports results even if they are lousy. That is sometimes helpful in finding better results.

@Axel Vogt Assuming that all parameters are positive we get
res1:=dsolve(Eq1) assuming positive;
series(rhs(res1),r=0,3); #So _C1 must be zero
diff(eval(%,_C1=0),r); #This is required to be zero

@mehdi jafari It is OK that you can opt out. But to have that as default seems strange to me.

Using the default:
restart;
with(CodeTools):
ode:=2*a*x^3*y(x)^3+2*x*y(x)+diff(y(x), x);
result:=Usage(dsolve(ode=0,y(x)));

you get the result from dsolve and a print of various info including cpu and real time. Do you need the time as actual output?

@mehdi jafari It used to be that if you answered a question or if you made a comment you were automatically signed up for email notification, i.e. subscription by default.

@samiyare Looking at the results obtained for NBT=1 (actually I think I looked at NBT=0.7) I saw that basically the graphs of phi, u, and T look like parabolas. So that is how I got those. You don't have to include the derivatives, as Maple can compute these, but I did initially since I was experimenting with some code of my own, which requires that the system be rewritten as a first order system. If the derivatives are included I assume those are used even if they are not really derivatives of the given approximations. For U and UF I used simple linear approximations.

About the case:

Phiavg:=0.1:
EPSILONE:=0.5:
Nr:=30:
lambda:=0.1;
Ha:=5;
NBT=0.2

I tried with those values, and with nothing else changed. I got no problem for NBT=0.2. In fact no problem from NBT=1 down to NBT=0.08. The first problem occurred at NBT=0.07.

I notice that the following device works:
latex(subs(erf=Erf,r));
and for the other example:
latex(subs(arctanh=Artanh,r));

The error comes up also with this simple version:
latex(erf(2*x));

but not with
latex(erf(x));
Same with arctanh.

The problem may be in `latex/print`, which is different from the procedure with the same name in Maple 17, where there is no such problem.
You may try
`latex/erf`(2*x);
p:=`latex/print`(erf(2*x)); #In Maple 17 p is a symbol, in Maple 18 a sequence
whattype~([p]); #Notice the type `*`
Compare with
p:=`latex/print`(sin(2*x)); #In Maple 17 and 18 p is a sequence
whattype~([p]);

I submitted an SCR about latex(erf(2*x)).

@Muhammad Ali The following two constructions both work. One uses piecewise, the other Heaviside:

pde:=diff(f(t,x),t)=diff(f(t,x),x,x)*piecewise(x>=1,1,0)+piecewise(x<1,1,0)*f(t,x);
#The alternative:
pde:=diff(f(t,x),t)=diff(f(t,x),x,x)*Heaviside(x-1)+Heaviside(1-x)*f(t,x);
res:=pdsolve(pde,{f(0,x)=piecewise(x>=1,(x-1)^2*x,0),f(t,0)=0,f(t,2)=2},numeric);
res:-plot3d(t=0..1);
p:=res:-value();
p(1,.5);
p(1,1.5);

Notice that f(t,x) does not appear inside piecewise or Heaviside.

@Muhammad Ali I'm not sure that I understand your example, because if f(t,x) is required to be 0 for x <= ln(b/K) then f(t,x) is known in that region, so that is that! In the remaning region x<ln(b/K) you must solve
diff(f(t,x),t)=diff(f(t,x),x,x) subject to initial and boundary conditions. The boundary condition at x=ln(b/K) ought to be chosen to be zero so that at least you have continuity.
For ln(b/K)=1 here is a simple example:
pde:=diff(f(t,x),t)=diff(f(t,x),x,x);
res:=pdsolve(pde,{f(0,x)=(x-1)^2*x,f(t,1)=0,f(t,2)=2},numeric);
res:-plot3d(t=0..1); p1:=%:
plots:-display(p1,plot3d(0,x=0..1,t=0..1));

@samiyare Notice that the results reported by LSSolve corresponds to the smallest values of a[1] and a[2], not to the last ones computed.

@samiyare I tried a different continuation but otherwise using the same loop as before:
N_bt:=cc*NBT+(1-cc)*1;#Works down to NBT=0.37
I also tried without continuation, i.e. removing the option continuation=cc and simply setting
N_bt:=NBT;
That worked down to NBT=0.41.

Your system is linear in x and y, so certainly you can trust solve to have found all solutions of this simple system. What does t have to do with this system?

@abscissa If you have several equilibrium points for a given set of concrete values of the parameters (i.e. in your case A, B, C, D, E, F, G, H, S, T, U, V, W, Y, x) then you may find that defining J as a function of the variables a, b, c, d, e, f, g, h, s, t, u, v may be convenient.
In two steps for clarity:
J1 := VectorCalculus:-Jacobian([eq_1, eq_2, eq_3, eq_4, eq_5, eq_6, eq_7, eq_8, eq_9, eq_10, eq_11, eq_12], [a, b, c, d, e, f, g, h, s, t, u, v]);
J:=unapply(J1,[a, b, c, d, e, f, g, h, s, t, u, v]):
#Now you can use J as a function as in
pt:=seq(1..12);
J(pt);