Unless you change the name of the imaginary unit the letter 'I' represents the imaginary unit.
By I*E did you actually just mean the name IE or EI?

Using Arrays (capital A) for both of Dpc and Z and continuing from the definition of 'vals':

#vals:={seq(seq(var[s][i]=s*i,i=1..NS),s=1..Ns)};
Z:=Array(1..Ns,1..NS,(s,i)->s*i); #Works as before with the changes below
## I also tried random numbers
r:=rand(0..99);
r();
Z:=Array(1..Ns,1..NS,(s,i)->r());
#Of course the elements of the Array could be defined one at a time as they are calculated.
#Anyway we got a 2-dim Array.
vals:={seq(seq(var[s][i]=Z[s,i],i=1..NS),s=1..Ns)};
#Proceeding as before:
eqs:=eval({seq(seq(seq(df_dz[s][i,j],i=1..NP),j=NP+1..NP+NC),s=1..Ns)},vals);
indets(eqs,function);
res:=solve(eqs,indets(eqs,function));
## Now defining Dpc to be a 3-dim Array. Default values will be zero
Dpc:=Array(1..Ns,1..NP,NP+1..NP+NC); #You could add the optional argument ,datatype=rational (if relevant).
res1:=evalindets(res,function,x->op(0,x)); #For convenience I remove the arguments
for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
Dpc[s,i,j]:=subs(res1,D[j-NP](z[s][i]));
od: od: od;

Dpc[1,1,3];
Dpc[1][1,3]; #Same as above

## Comment The only arrays or Arrays appearing are Dpc and Z. Quantities such as phi, phi[1] etc. are implicitly defined as tables by the assignments like phi[1][2]:=..... In fact phi is a table of table(s) (plural if Ns>1).
eval(phi);

@tira Was my guess that you were missing a parenthesis in f correct?

Specifically was my guess that f should have been

f:=((1-exp((-t*u^2)/(1+alpha*u^2)))*sin(y*u))/(u^3);

and NOT

f:=(1-exp((-t*u^2)/(1+alpha*u^2))*sin(y*u))/(u^3);

correct or not?

@Carl Love The coeffs option is used in dsolve with type=formal_series (and in type=formal_solution).

Let the ode be called 'ode'. We consider the homogeneous equation.

dsolve(eval(ode,q=0),w(x),'formal_series');
res:=dsolve(eval(ode,q=0),w(x),'formal_series','coeffs'='mhypergeom'); #The most interesting (see below)
dsolve(eval(ode,q=0),w(x),'formal_series','coeffs'='hypergeom');
dsolve(eval(ode,q=0),w(x),'formal_series','coeffs'='polynomial'); #NULL
dsolve(eval(ode,q=0),w(x),'formal_series','coeffs'='rational'); #NULL

plot(eval(rhs~([res]),{_C1=0,_C2=0,_C3=1,Nr=1,Nl=2,EJ=1}),x=0..15);

In the following code I avoided using alias. I didn't need to use unapply after all.

restart;
printlevel:=3:
Ns:=1:
NP:=2:
NC:=2:
NS:=NP+NC:
for s from 1 to Ns do
  indepvar[s]:=seq(z[s][j],j=NP+1..NP+NC);
  depvar[s]:=seq(z[s][j](indepvar[s]),j=1..NP);
  var[s]:=depvar[s],indepvar[s]
end do;

for s from 1 to Ns do
for i from 1 to NP do
phi[s][i]:=add(var[s][j]^2,j=1..NS)
od: od;

for s from 1 to Ns do
for i from 1 to NP do
f[s][i]:=expand(depvar[s][i]-phi[s][i]);
od: od;

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=diff(f[s][i],z[s][j]);
df_dz[s][i,j]:=convert(df_dz[s][i,j],D)
od: od: od;

vals:={seq(seq(var[s][i]=s*i,i=1..NS),s=1..Ns)};
eqs:=eval({seq(seq(seq(df_dz[s][i,j],i=1..NP),j=NP+1..NP+NC),s=1..Ns)},vals);
indets(eqs,function);
solve(eqs,indets(eqs,function));

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