@leoteo Let me say right away that I don't see how we could find a function h(x,y) which satisfies both EQ1 and EQ2 no matter what initial or boundary conditions you impose. I see one trivial exception to this: h(x,y)=L:

eval([EQ1,EQ2],h(x,y)=k);

So if k=L both equations are satisfied. Initial and boundary conditions would have to conform to this of course.

Below I only consider EQ1, but be aware that I determine h(x,y) completely from EQ1 and the initial conditions I choose below. Thus either EQ2 is satisfied for that solution or it is not.

Let me give you the whole code. Notice though that I haven't corrected the squaring errors. You can just do that. But try this first.

restart;
eq1:=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;
EQ1:=subs(U2(x,y)=diff(U1(x, y), y),U1=h,eq1);
eq2:=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;
EQ2:=subs(U3(x,y)=diff(U1(x, y), x),U2(x,y)=diff(U1(x, y), y),U1=h,eq2);
#I'm not using your conditions as they make no sense to me as I mentioned:
#conds:={h(x0,y0)=K1,D[1](h)(x0,y0)=K2,D[2](h)(x0,y0)=K3};
#But I'm using these:
#Example initial conditions.
ic1:={h(x,0)=sin(x),D[2](h)(x,0)=1};
#You would hope that this would work:
pdsolve(eval(EQ1,{L=1,D2=1,alpha=1}),ic1,time=y,numeric,range=0..Pi);
#but error.
#Boundary conditions shouldn't be necessary since the initial conditions clearly determine u(x,y). #However, picking h(0,y)=y is consistent with the initial conditions D[2](h)(x,0)=1 and h(x,0)= sin(x).
pdsolve(eval(EQ1,{L=1,D2=1,alpha=1}),ic1 union {h(0,y)=y},time=y,numeric,range=0..Pi);
#No luck, and this error message seems to contradict the other. A royal round around.
#The pde EQ1 is probably too far from a standard form:
#Quote from help page for pdsolve,numeric:
# "The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution."
#The problem here is that it is too simple: it really is just an ode in y for each fixed x. I would consider this failure a weakness or a minor bug.
##SOLUTION 1. Tricking pdsolve into believing that the pde is more of a standard pde.
#We add a harmless term involving a derivative w.r.t. x. to make pdsolve see the pde as more of a ##standard one. The term will be totally harmless because it is zero in our chosen interval x=0..Pi, in #fact it is zero for x < 5.
EQ1a:=EQ1+(0=piecewise(x<5,0,1)*diff(h(x,y),x));
res:=pdsolve(eval(EQ1a,{L=1,D2=1,alpha=1}),ic1 union {h(0,y)=y},time=y,numeric,range=0..evalf(Pi));
res:-plot3d(y=0..1);
###########
##SOLUTION 2. Treating EQ1 as an ordinary differential equation in y for each fixed x. I shall be using the parameters option in dsolve. The parameter being x.
EQ1d:=subs(h(x,y)=h(y),EQ1);
res:=dsolve({eval(EQ1d,{L=1,D2=1,alpha=1}),h(0)=sin(x),D(h(0)=1},
numeric,output=listprocedure,parameters=[x]);
H,H1:=op(subs(res,[h(y),diff(h(y),y)])); #Extracting procedures for h(y) and h'(y).
p:=proc(x,y) res(parameters=[x]); H(y) end proc;
p(0.123,1); #Checking p on (x,y)=(0.123, 1).
plot3d(p,0..Pi,0..1);

@leoteo The corrections don't change the basic problem: You have more equations than dependent variables and the initial/boundary conditions are given at just one point.

In the original h-formulation you have 2 equations but one dependent variable (h).
In the systems version you have 4 equations and 3 dependent variables (U1,U2,U3).

There is no real reason for rewriting the equations. If Maple could handle the problem it would do the rewriting itself.

Incidentally, you still have a discrepancy between the systems version and the original: The cosines are squared in the systems version, not in the h-version. Correcting this won't affect the basic problem.
Please notice that I edited my answer above.

@leoleo000 I notice right away a few issues in your procedure errorf:

Periods instead of commas in the local declaration.
A missing parenthesis in the definition of ef.
x is declared local, yet it is used as a formal parameter.
y and z are also formal parameters but are not used in the procedure.
a,b,c are used in the procedure. What are they?

Don't you want the formal parameters besides f to be a,b,c?

So what is the problem and why is it so urgent?

@emmantop But bd represents the fourth derivative, not the second derivative, so it cannot be used in eq[12].

To illustrate that we can do:
restart;
bd := (f[i]-4*f[i-1]+6*f[i-2]-4*f[i-3]+f[i-4])/h^4;
bd0:=eval(bd,i=2);
for i from -2 to 2 do f[i]:=(i*h)^4 end do; # Corresponding to f(x)=x^4
bd0; #Gives the correct fourth derivative of x^4.
for i from -2 to 2 do f[i]:=(i*h)^2 end do; # Corresponding to f(x)=x^2
(f[1]-2*f[0]+f[-1])/h^2; #This gives the correct second derivative of x^2
bd0*h^2; # This does not

@smith_alpha You could add the .txt by using  cat("D:/",file_name,".txt").

If you really defined f as a procedure, then plot(f,-2..2) should work.

What was your definition of f (exactly)?

@John Starrett With Digits:=30 as I suggested earlier I cannot see any difference between 
display(Array([PW1a,PW2a,PW3a])) and display(Array([PW1,PW2,PW3])).

Before computing DEP1 and DEP1a I set Digits back to 10 and used maxfun=0 in both.
Again no difference.

Note:  I used the first version of your posted worksheet.

@emmantop I would linearize after discretization. By linearization I simply mean using Newton's method on the system of equations ensuing from discretization. That involves a starting guess and iteration like in the one dimensional Newton's method.
However, I would begin by writing the ode-system as a first order system. Then discretize.

Obviously you would have to replace infinity by some number like 10 and use a midpoint method.
To test your results you could compare with

res:=dsolve({ode1,ode2} union subs(infinity=10,{bc1,bc2}),numeric,method=bvp[midrich]);

@adel-00 
The straightforward approach actually works, you just have to interpret the result correctly (in the first version of this comment I didn't have it straight).
eqs:=subs(x(t)=x,y(t)=y,z(t)=z,rhs~(dsys))=~0;
res:=solve(eqs,{x,y,z});
#x=RootOf(conjugate(_Z)+_Z) means that x is a number whose real part is 0.
#Thus x = I*x2, where x2 is real.
#The result is therefore
evalc(subs(RootOf(conjugate(_Z)+_Z)=I*x2,res));
#i.e. there are infinitely many solutions parametrized by x2.
#We check these solutions below.
#We could split instead:
eqs2:=subs(x(t)=x1+I*x2,y(t)=y1+I*y2,z(t)=z,rhs~(dsys))=~0;
eqs2R:=(evalc@Re)~(eqs2);
eqs2I:=(evalc@Im)~(eqs2);
solve(eqs2R union eqs2I,{x1,x2,y1,y2,z});
#Notice that the solutions are the same as found directly above.
#Checking:
eval(eqs2,%);
simplify(%) assuming real;
 

Please note that I revised my second answer totally. I submitted a bug report expressing this conflict about plotting before saving.

@adel-00 I think you meant 0.25*z(t)^2+Re(y(t))^2+Im(y(t))^2 because that quantity should be and appears to be (roughly) constant (=1/16=0.0625), whereas the one you actually have doesn't:

plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2+Im(y(t))^2 ],0..10);
plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2+Im(y(t)) ],0..10);

By saying "I have written a demonstration worksheet to show this problem", I presume that you intended to include that worksheet in your question, but I don't see any link to it.

@adel-00 If you just want to plot the quantity 0.25*z(t)^2+Re(y(t))^2 then you can just use res from before:

plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2],0..2);
# but if you want to get actual numbers you can use output=listprocedure, which would be OK to use in any case.
#Notice that this time I have closed the assignments to resL, Y, Z by a colon. If you replace colon with semicolon you will know why (a huge output: it is usually hidden. It is surely not intended. I shall report it as a minor bug).

resL:=dsolve(dsys union {x(0)=0,y(0)=0,z(0)=-1/2},numeric,output=listprocedure):
#You can use odeplot as before:
plots:-odeplot(resL,[t,0.25*z(t)^2+Re(y(t))^2],0..2);
Y,Z:=op(subs(resL,[y(t),z(t)])):
Y(0.12345);
Z(0.12345);
#You could if you wish define the quantity as W:
W:=t->0.25*Z(t)^2+Re(Y(t))^2;
W(.987654);
#And you can plot W (same plot as above):
plot(W(t),t=0..2);

@RAfossey Here is a stepwise solution by laplace transform. As you are doing I'm assigning to a,b,c,d immediately although that doesn't seem to be what the problem formulation called for.

restart;

with(inttrans):
a:=1;b:=1;c:=1;d:=1;
sys:=diff(x(t),t,t)=a*x(t)+b*y(t),diff(y(t),t,t)=c*x(t)+d*y(t);
ics:=x(0)=1,D(x)(0)=0,y(0)=0,D(y)(0)=1;
#Applying laplace elementwise (notice the tilde ~ ) to the system:
laplace~([sys],t,s);
eval(%,{ics});
#Solving for the laplace transforms:
L:=solve(%,{laplace(x(t),t,s),laplace(y(t),t,s)});
#Elementwise inverse transformation
sol:=invlaplace~(L,s,t);
#Extracting the solutions:
X,Y:=op(subs(sol,[x(t),y(t)]));
#Plots:
plot([X,Y],t=0..2);
#This is a parametric plot:
plot([X,Y,t=0..2]);