@emmantop When looking into the errors by using the result from dsolve as the correct result (not unreasonable) it is striking how bad they are. They indicate that the overall method is first order in h, i.e. a halving of the stepsize only gives a halving of the error.
The culprit is the discrete version of the first derivative used at the boundary:
eq[N] := (f[1]-f[0])/h = 1
and
eq[N+1] := (f[N]-f[N-1])/h = 0

If you replace those by by the symmetrical differences:
eq[N] := (f[1]-f[-1])/h/2 = 1;
eq[N+1] := (f[N+1]-f[N-1])/h/2 = 0;

you get a second order method.

I have uploaded a new version of the worksheet. It includes additional improvements, in particular to the loops, where fsolve is now taking advantage of previously found results in the guess. 

MaplePrimes15-02-04odebvp_discrete_Update.mw

The loops in the worksheet now goes through N= 5,10,20,40,80.
A time to stop would be dictated by the tolerance given. If that is 2e-4 then the actual error (as compared to dsolve) is less than that for N=80, but not for N=40.
Without using a "correct" result you can find the differences of results for 5 and 10, 10 and 20, 20 and 40, 40 and 80. When a difference gets below the tolerance you could stop. If the tolerance is 2e-4 as above you would have go on to N=160.
To improve results already found for e.g. N=40 and N=80 you could use Richardson extrapolation knowing that the order of the method is 2.

@Carl Love You are right. I removed the exception in my answer above.

In my earliest version I had
int(sin(y*u)/u^3, u = 0 .. 1);

in which case y=0 is an exception.

@wo0olf That pde is quite different from any interpretation of the "2" in my comment above! Now at least you are getting a solution, so that is not the problem.

@wo0olf Does

diff(T(x, y), x)2*(diff(T(x, y), y, y))

mean

diff(T(x, y), x)^2*(diff(T(x, y), y, y))

or

diff(T(x, y), x)*2*(diff(T(x, y), y, y))

or shouldn't the number 2 be there at all?

I cannot add anything to the description in the help pages for pdsolve/numeric. Wanting a complete description of the method being used is not encouraging people to contribute with anything.

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