@Preben Alsholm  

Thanks for your quick reply!!! Your code is running but the final result shows only the procedure: Using rkf45 to solve the first ODE. I have two questions: 

1. Could I show the final result function of h(x,y) or h'(x,y)??

2. How about the second PDE? I just kept the same condition from the first ODE and solve second PDE but something is still not right... see the below code:

res := pdsolve({eval(EQ2, {D2 = 1, L = 1, alpha = 1}), h(0) = 0, (D(h))(0) = 1}, numeric);
Error, (in pdsolve/numeric) invalid subscript selector

I am just beginner, sorry to ask you too many simple questions for you.

@Preben Alsholm Hello, 

Sorry, because I forgot to update the new model of my PDEs. The two PDEs could be rewritten as:

https://drive.google.com/file/d/0B141s-K48fmVR3p4SGFhWjR4VFE/view?usp=sharing

P/S: the last conmmand was wrong when I copied it to this website. The right one is:

pdsolve(PDEs, IBC, (U1, U2, U3), nummeric); 

with empty response.

I already changed the input type and add more initial condition and used pdsolve () function but the response is empty:

My code:

with(PDEtools); declare(U1(x, y), U2(x, y), U3(x, y));

PDEs := [diff(U1(x, y), y) = U2(x, y),

diff(U1(x, y), x) = U3(x, y), (diff(U2(x, y), y))*((L-U1(x, y))^2+(D2+tan(alpha)*(L+U1(x, y)))^2)+.5*U2(x, y)*(2*tan(alpha)*U1(x, y)-D2) = 0,

(diff(U3(x, y), x))*((L-U1(x, y))^2+D1^2*cos(alpha+2*U2(x, y))^2)+.5*U3(x, y)*D1*cos(alpha+2*U2(x, y))+.5*D1*sin(alpha+2*U2(x, y))*(L-U1(x, y)) = 0];

IBC = {U1(0, 0) = 0, U2(0, 0) = 0, U3(0, 0) = 0};

pdsolve(PDEs, IBC, U1, U2, U3, nummeric);

What should I do next???