@Chia Infinitely many solutions are returned by dsolve(de); For each constant _C1 you got a solution.
If you add an initial condition like y(0)=35, then there is exactly one solution (under very general conditions on f). Thus if the operator T is defined on (say) C^1(R) (functions differentiable on R with continuous derivative) by T(y) = y' and satisfying y(0)=35, then the equation T(y)=f has a unique solution y for given f as long as f is continuous on R. Thus the inverse of T, T^(-1), exists and is defined on C(R).
The main point to be made, though, is that this doesn't help us one bit in finding the solution.

The obvious first answer is that D[1, 1]+D[1, 2]+D[2, 2] doesn't have an inverse.
Just think of the equation y'(x)=f(x), where y is the unknown. Using D you have D(y)=f. Clearly, if there is a solution,there are infinitely many. To single out one of those you need an extra condition. The same holds for the differential operator D[1, 1]+D[1, 2]+D[2, 2].

@wo0olf But D[1](M)(x,0)=1000 when diff(T(x,y),y)=M(x,y) means the same as D[1,2](T)(x,0)=1000, i.e. the mixed derivative of second order.
Notice that in Maple diff(T(x,y),y) = D[2](T)(x,y).

@lijiwei You must correct (or reformulate) the problem, in particular the definition of Y.
Is h(x) always supposed to nonnegative? If h(x)>0 then s(x)>1, thus (1-s(x)) < 0 so that (1-s(x))^(-1+1/n) = 1/(1-s(x))^(2/7) is imaginary.

@wo0olf Well, then you got the problem that the boundary conditions can only contain derivatives normal to the boundary, thus D[1](T)(x,0)=1000 and D[1](T)(x,0.000025)=2000 are not acceptable.

Could you provide text instead of an image, please.

You have
Y := x->f-(1/2-(1/2)/n)*ln(s(x))+2*ln(1-(1-s(x))^(-1+1/n));
and with the values of your parameters then that evaluates to

thus you have a problem even though you only need diff(Y(x),x):
At the one end x=0 you want h(0)=0. If h(x)>0 for some x then you get an imaginary result because of the last term in Y, try
eval(Y(x),{diff(h(x),x)=1,h(x)=-.1});
If on the other hand h(x)<0 for some x then the fact that you have h(x)^n gives an imaginary result, try
eval(h(x)^n,h(x)=-.1);

The error I get from dsolve is not the one you report, but:
Error, (in fproc) unable to store 'HFloat(HFloat(undefined))+HFloat(HFloat(undefined))*I' when datatype=float[8]

@wo0olf BC1, BC2, and IC already determine a solution:

restart;
pde:= diff(T(x, y), x)-1.775*10^(-5)*(diff(T(x, y), y, y))/(-1.218493440*10^11*y^2+4.244913600*10^6*y+0.33e-1) = 0;
sol:=pdsolve(pde,{T(0, y) = 0,D[2](T)(x, 0) = 1000, D[2](T)(x,0.000025) = 2000},numeric);
sol:-animate(x=1e-3,frames=100);
sol:-plot3d(x=0..1e-3);

@wo0olf This is a new problem: The pde is different and the conditions are different.
You cannot impose conditions like (D[1](T))(0, y) = 1000, (D[1](T))(x, 0) = 2000. The first one because the pde is first order in x, the second because the boundary conditions (as stated in the error message) can only contain derivatives normal to the boundary, thus (D[2](T))(x, 0) = 2000 would have been acceptable.
An example of a problem that can be handled:
sol:=pdsolve(pde,{T(0, y) = 0,D[2](T)(x, 0) = 1000, D[2](T)(x,1) = 2000},numeric);
Maybe I should add that D[1](T) means the derivative of T w.r.t. the first variable (called x in your case). Similarly D[2](T) means the derivative of T w.r.t. the second variable, i.e. y.

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