@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

@necron Now with names (and more) changed I again suggest that before attempting to plot you should try to see if you can get values out of quantities like

evalf(g2(2.345)); No chance if you cannot do these:
evalf(eval(c2,{u=1,r=2})); #Comes up with a number
evalf(eval(c2,{u=1,r=20})); #This doesn't
evalf(eval(c3,{u=1,r=2})); #This doesn't either

@Prashanth I think this little session shows what happens. I'm using the same example as I did above.

restart; #Important
p := proc()
    with( ListTools );
    Reverse( [ 1, 2, 3 ] )
end proc:
eval(Reverse); #with(ListTools) has not been executed, but Reverse inside is now the global :-Reverse
p(); #This will cause execution of with(ListTools)

eval(Reverse); #This is ListTools:-Reverse
Reverse([1,2,3]); #It works
p(); #Reverse inside p remains :-Reverse
op(0,p()); # This is :-Reverse
%-Reverse; # Not zero
###Now executing without a restart the same procedure definition once more
p := proc()
    with( ListTools );
    Reverse( [ 1, 2, 3 ] )
end proc:
eval(Reverse); #with(ListTools) was already executed above so the Reverse seen is ListTools:-Reverse
p();
###If you try the whole thing from the start and use another name (q) for the procedure the second time, then you will have both p and q. p still doesn't work, but q does.

You are missing a value for 'u' for one thing.

Since you have an iterated integral you ought to see if you can get a value out of just e.g. h2(1.2345) with the value of u that you have. If that fails then forget about plotting.

@Marcel123 My general advice is: Do not ever assign to the variables used in odes, in your case x and t.
The variables are used in the solution as just variables and shouldn't be changed.

Use other letters, as in

X:=t->rhs(sol(t)[2]);   ## Capital X

@Z1493 From what I have seen so far my answer to the question "I only wanted to know if I have a PDE system with composite functions (as in the above case, f=f(x(t),t) ) how would I proceed to integrate it"
would be: Solve the pde system, then insert x=x(t).
But I think that what is needed now is an actual physical situation where this comes up.

@Z1493 Let me point out that
(1) The pdes are not the same.
(2) The solutions given by pdsolve are building blocks. See my previous reply with the infinite sum.

Question: Are you asking if given one system of (linear) pdes with coefficients depending on t and another where t is replaced by ln(x) in the cofficients only then the general solution for the two different systems will agree if we replace x by exp(t) in both?
Or expressed differently: Along the curve (x,t) = (exp(t),t) the solutions are the same?
That may be right, but is it interesting? Because of my point (2) above this is not disproved by your example (where you don't really test for that though: You need to replace x by exp(t)).

@Z1493 It is not clear to me what it is you want to solve.
So instead of your abstract formulation with A and B could you formulate the problem (i.e. the equations) concretely from the start.
You don't mean (I take it):

restart;
eq1:=D[2](f1)(exp(t), t) = f2(exp(t), t);
eq2:=D[1](f2)(exp(t), t) = -f1(exp(t), t);

and then hope to solve this for f1 and f2?

@Z1493 Your pde system could be rewritten as a second order pde in f1 (which I shall call f):
pde:=diff(f(x,t),x,t)+f(x,t)=0;
pdsolve(pde);
pdsolve(pde,build);
combine(%);
#By linearity any (convergent) infinite linear combination of solutions is again a solution:
S:=Sum(c(n)*exp(q(n)*x-t/q(n)),n=1..infinity);
pdetest(f(x,t)=S,pde); #Gives 0 as expected.
##Example
eval(S,{c(n)=1/n^4,q(n)=I*n}); #Picked c(n) for fast convergence
solR,solI:=op(evalc([Re,Im](%))); #If you have a complex solution then the real and imaginary parts are also solutions:
pdetest(f(x,t)=solR,pde); #Similar for solI
#For plotting purposes we can replace infinity by e.g. 10:
plots:-animate(plot,[eval([solR,solI],infinity=10),x=0..10],t=0..10,frames=100,size=[1000,400]);

We must expect that there are other nontrivial solutions corresponding to other values of L.
To find one more I found it easier to use a 2nd order ode satisfied by y(xi).

UU:=(isolate(sys[2],U(xi)));
eval(sys[1],UU);
ode:=collect((lhs-rhs)(%),diff,factor)=0;
odeL:=subs(y(xi)=L*y(xi),convert(ode,D));
resL:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=1,y(xi)=(1-xi)*(xi-0.8)]);
resL(.9); #Finds L = 1 approximately (as before)
UU; #(the right hand side can be plotted)
plots:-odeplot(resL,[[xi,y(xi)],[xi,rhs(UU)]]);
## Now looking for another value of L near 4:
resL2:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=4,y(xi)=(1-xi)*(xi-0.8)*(xi-0.9)]);
resL2(1); #L very close to 4
plots:-odeplot(resL2,[[xi,y(xi)],[xi,rhs(UU)]]);

@RomchikM In the code I gave in my answer above you will see the line
plots:-odeplot(res,[[xi,y(xi)],[xi,U(xi)]]);

This shows the graphs of y and U.

########
The problem you are facing is very similar to the following problem.

restart;
ode:=diff(y(x),x,x)+one*y(x)=0;
bcs:= y(0)=0, y(Pi)=0;

dsolve({eval(ode,one=1),bcs}); #infinitely many solutions (as is well known)
dsolve({eval(ode,one=1+1*10^(-8)),bcs}); #No nontrivial solution
res:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric); #Numerical attempt
plots:-odeplot(res,thickness=3);# same result 0
res2:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric,approxsoln=[y(x)=sin(x)]);
plots:-odeplot(res2,thickness=3); #Tiny, but it resembles a*sin(x) for a very small a
resL:=dsolve({ode,bcs,D(y)(0)=1},numeric);
resL(0); #'one' is very close to 1
plots:-odeplot(resL);



So I would suggest using the results I gave in my answer with the adjusted parameter L. Your input equations contain floats and are surely rounded, so why not?

@MTata94 Sure.
But you could also define the procedure q outside p, but then it shouldn't be declared local to p.
If you have a lot of subroutines used in p you may want to define them outside p and then wrap the whole thing in a module where the subroutines could be declared local. That is another story.