@Axel Vogt Epsilon prints like a capital epsilon, which ought to be fine.

By the way I don't see the solution found by Carl in your reduced system.

Another thing: I cannot make fsolve give a result on Sys[2] without omitting convert/rational. Without that, no problem.

Using solve on the equations without assignment to the constants gives two real solutions altogether when assigning afterwards:

For ppsi = .1561816797:
{A1 = .2178679378, A2 = 0, C0 = 0, C1 = 0, C2 = 0, E0 = 2.923409461, H = 1.169518511, H1 = -3.370336272, K = .2120924177, L0 = -1.923409459, L1 = 1, R = 0.1077611106e-1, W = .4031652600, Y = .6858871836, Epsilon1 = -1.108403687}
and
{A1 = .1151377592, A2 = -0.2771025818e-1, C0 = -0.2106664520e-1, C1 = -0.3108676424e-1, C2 = -0.4587284305e-1, E0 = 2.923409461, H = .9700004405, H1 = -.4209490949, K = 0.8582496732e-1, L0 = -3.093516599, L1 = -.1854153076, R = .655446246, W = .3221564949, Y = .4545705612, Epsilon1 = .2686838282}

For ppsi=1.947717707:
{A1 = 0.1146472975e-1, A2 = 0, C0 = 0, C1 = 0, C2 = 0, E0 = .2344194943, H = 0.9509540701e-1, H1 = 0.6343537741e-1, K = 0.1116080812e-1, L0 = .7655805057, L1 = 1, R = .3632312348, W = .3518971259, Y = 0.4867856263e-1, Epsilon1 = 0.3694572998e-1}
and
{A1 = 0.1429702382e-2, A2 = 0.5359982543e-1, C0 = 0.4074911503e-1, C1 = 0.6013098524e-1, C2 = 0.8873162983e-1, E0 = .2344194943, H = .5898272830, H1 = 1.456583019, K = 0.5218750955e-1, L0 = .1114825845, L1 = .3373447074, R = .655446246, W = .3221564946, Y = .2764103064, Epsilon1 = 0.4295105618e-1}

The above results were computed with Digits=10, but confirmed with Digits=15.

(NOTE: Second value of ppsi).
Since you know the solution you can use that as a starting point for fsolve just to verify the result, always a good idea.
With the method I mentioned in my previous comment (but didn't elaborate on) I found

{A1 = 0.142970229228308e-2, A2 = 0.535998252919882e-1, C0 = 0.407491149283945e-1, C1 = 0.601309851312542e-1, C2 = 0.887316296122945e-1, E0 = .234419494359216, H = .589827281889861, H1 = 1.45658301648039, K = 0.521875094397844e-1, L0 = .111482584275646, L1 = .337344707782075, R = .655446246865986, W = .322156494043605, Y = .276410305860281, Epsilon1 = 0.429510561240001e-1}

which is the one you are talking about.

Using a starting point can help fsolve quite a lot. So any method that will provide you with an approximate solution may help.

I was using your system as a test on some solving method I'm looking into and found one more solution (there may be more):
(NOTE: First value of ppsi).
It has  
sol1:={A2=0,C0=0,C1=0,C2=0,L1=1};
#and
sol2:={A1 = .217867938869731, E0 = 2.92340946074119, H = 1.16951851156081, H1 = -3.37033625373258, K = .212092418616612, L0 = -1.92340946074119, R = 0.107761071012707e-1, W = .403165260852399, Y = .685887185020117, Epsilon1 = -1.10840368196766};

#You can test the soluion with fsolve by doing

eval(eq,{A2=0,C0=0,C1=0,C2=0,L1=1});
neweq:=select(hastype,%,name);
fsolve(neweq,sol2);
#However, the following gave me the solution found by Carl (I'm also using Digits=15):

fsolve(eq, sol1 union sol2);

@9009134 I'm looking at your latest worksheet called frekans.mw. I don't know what you call the original.

1. I don't see any code like that in that worksheet.
2. The try ... end try construction is used when you want to try something, but expect that some error might occur and don't want that to stop the process from running. Here you want to try different extra bcs.
If dsolve doesn't work with that extra bcs (b) and results in an error then you want that error caught (the catch clause). In case of error the last error is printed, but the loop is not interrupted, you go on to the next b in extra_bcs.
3. Indeed! Try changing it.
4. b is a loop variable just as in this simple construction:
   for b in [7,9,13] do b^2 end do;
4 II. Only the b's that didn't result in errors produce results, so you can pick any of those.
5. Since RES was defined as res[b] for one of the successful extra bcs b, RES(0) will give you the value of the solution for theta=0. That was just intended for getting the value of omega, which supposedly is constant. For that worksheet I get 9.03506039191669.
You might as well have done RES(1);

@9009134 dsys4 uses theta. You use q later. Change each q to theta (or vice versa).

@Kitonum SCR stands for "Software Change Request" mostly thought of as a euphemism for bug report.
You can find a way to do it in the menu line on top on this site under "More".

You start with a substitution into Ca[2] and receive an answer involving sin. Where did that come from?

If we are to help you we must see the code from the start. You can forget about the output: remove it before copying.

@Al86 Just do:

eval~(Ymono,t=~T);

@rit

I don't know what you mean by "I choose one or zero". One or zero where?

The RootOf result from solve will give 3 different sets of solutions for (a,b,c). Choose the one that gives you real solutions in a neignorhood of t=1 (if real solutions are what you want).
The other two solutions will be imaginary.
Forget about symbolic solutions.
Example:

eval[recurse](convert({ode},D),{t=1,D(b)(1)=0}); #Example
solve(%,{a(1),b(1),c(1)});
## So it looks like we can get away with:
ics1:=a(1)=1,b(1)=0,c(1)=1; #Example
res:=solve({ode},diff({a,b,c}(t),t));
Ares:=allvalues(res):
nops([Ares]);
Ares[1];
sol1:=dsolve(Ares[1] union {a(1)=1,b(1)=1,c(1)=1},numeric);
plots:-odeplot(sol1,[[t,a(t)],[t,b(t)],[t,c(t)]],0.5..1.62);
plots:-odeplot(sol1,[t,a(t)],0.5..1.62,thickness=2); p1:=%:
sol2:=dsolve(Ares[2] union {a(1)=1,b(1)=1,c(1)=1},numeric);
sol3:=dsolve(Ares[3] union {a(1)=1,b(1)=1,c(1)=1},numeric);
sol2(1.2345);
sol3(1.2345);
plots:-odeplot(sol2,[[t,Re(a(t))],[t,Im(a(t))]],-1..5,color=[green,red]);p2:=%:
plots:-odeplot(sol3,[[t,Re(a(t))],[t,Im(a(t))]],-1..5,color=[green,red]); p3:=%:
plots:-display(p1,p2,p3);

@taro The level of evaluation inside a procedure has nothing to do with the fact that your procedure g:=x->f1 doesn't work.
The reason is that before any evaluation at all takes place (no matter the level) when you do g(2); the variable x used in x->....  is replaced by 2 (simple substitution). Now no x is seen in the body. It will be seen after evaluation surely, but that will take place after the substitution and that is too late, so that x doesn't get the value 2.

I tried using the first part of your code till you have defined all the equations and the variables, eqs and vars, respectively.

Then I used fsolve:

res_phi:=fsolve(eqs,{vars[]}=~1); #Giving start values 1 to all variables

After that I defined the matrix M by

M:=Matrix(subs(res_phi,[seq([seq(phi[p,q],p=0..Numy+1)],q=0..Numz+1)]));

This matrix I compared to your matrix sol3. The difference matrix M-sol3 had elements of absolute value less than 5e-11:

max(abs~(M-sol3));
     4.98876495669264841*10^(-11)

@Carl Love As for second order in time we can take the wave equation, which of course is not elliptic, but you said "any PDE that isn't first-order in time".

We can compare with the exact solution. Incidentally the exact solution found by pdsolve is wrong!

restart;
pde:=diff(u(x,t),t,t)=diff(u(x,t),x,x);
ibcs:=u(x,0)=0,D[2](u)(x,0)=sin(x),u(0,t)=0,u(Pi,t)=0;
sol:=pdsolve({pde,ibcs}); #WRONG! But we can see that u(x,t) = sin(x)*sin(t) is a solution.
pdetest(u(x,t)=sin(x)*sin(t),[pde,ibcs]); #Check
eval((rhs-lhs)~({pde,ibcs}),u=((x,t)->sin(x)*sin(t))); #Double check
res:=pdsolve(pde,{ibcs},numeric); #Numerical solution
res:-animate(t=4*Pi,frames=100); p1:=%:
plots:-animate(plot,[sin(x)*sin(t),x=0..Pi],t=0..4*Pi,frames=100); p2:=%:
plots:-display(p1,p2); #Doesn't look bad at all!

@Carl Love I try, but often fail, to choose my wording carefully.
I wrote
"Thus the IBC given by
IBC2:={u(0, y) = 1, D[1](u)(1, y) = 0, u(x, 0) = 1, D[2](u)(x, 0) = 0};
would not be rejected by pdsolve/numeric." (emphasis added).
That was what I meant.

##############
It may be fair to say that the time-based solver is not doing a great job on elliptic pdes, but it isn't failing entirely.

Take the laplace equation:

restart;
pde := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0;
U:=evalc(Re(exp(x+I*y))); #Harmonic, so solves pde
#u(x,y) = U satisfies
ibcs:={u(0, y) = cos(y), D[1](u)(1, y) = exp(1)*cos(y), u(x, 0) = exp(x), D[2](u)(x, 0) = 0};
pds:= pdsolve(pde,ibcs, numeric,timestep=0.001);
pds:-plot3d(y=0..0.3); p1:=%:  #Yes, don't go too far with y. Already a ripple starting.
plot3d(U,x=0..1,y=0..0.3); p2:=%:
plots:-display(p1,p2);
pds:-value()(.2345,0.07913);
eval(U,{x=.2345,y=0.07913});

Let me add that the error message may be misleading. It is not so much a question whether the pde itself is elliptic or not. You have to include the boundary conditions, IBC in your case. They are given on the square [0,1]x[0,1]. Thus the problem is not handled by the algorithm used by pdsolve/numeric.
This algorithm is time based, which basically means that it starts from time t=t0 with enough information to march forward step by step. In your case time could be taken as y. Since the pde is second order in time (y) in order to get going the values of u(x,y0) and D[2](u)(x,y0) should be known (for some y0,here 0, and all x in some interval, here [0,1]).

Thus the IBC given by
IBC2:={u(0, y) = 1, D[1](u)(1, y) = 0, u(x, 0) = 1, D[2](u)(x, 0) = 0};
would not be rejected by pdsolve/numeric. You can give the optional argument 'time'=y, but it is not necessary since it can be deduced from IBC2.

First you have to correct a few syntax errors: missing arguments in 'u' in two places, a space between D[i](u) and what follows in the IBC (disastrous in 2D-math input: avoid it like the plague!).

After that you got:

PDE := diff(u(x, y), x, x)+diff(u(x, y), y, y) = u(x, y)^(1/2)+diff(u(x, y), x)^2/u(x, y)^(3/2);
IBC := {u(0, y) = 1, u(x, 0) = 1, D[1](u)(1, y) = 0, D[2](u)(x, 1) = 0};
pds := pdsolve(PDE, IBC, type = numeric);

The result is the error message:

Error, (in pdsolve/numeric) unable to handle elliptic PDEs