@shabahang I believe you are wrong. Although pi and Pi look the same in output when entered as
pi;
Pi;

they are not the same. pi is just a name, Pi is the mathematical constant with approximate value 3.14159265.

Try
evalf(pi);
evalf(Pi);

Could you give us the lines of code you are using?

@ger89 I shall try to answer the questions one by one.

1. You don't know right away, but by trying first without the change _c[1]=c^4*L^(-4), as you would surely do, you may see that that change would make results look much nicer and more compact.

2. As I wrote earlier surely you didn't mean I*E, but EI or IE, i.e. a name. I was once exposed to some engineering and vaguely remember seeing EI used in this context (so I should have used EI, not IE).
Anyway IE or EI is just a factor in what you wrote as
eval(IE*(D[1, 1](w))(x, t), x = L) = 0 (I have replaced I*E by IE)
and presumably IE<>0 so that is equivalent to
eval(D[1, 1](w)(x, t), x = L) = 0
which is the same as
D[1, 1](w))(L, t) = 0

3. My answer here is similar to my answer to question 1. You observe that besides parameter c the equation resulting from
eval(ode1,{_c[1]=c^4*L^(-4),_F2=G});
includes IE/L^4/m. Thus it only depends on IE/m, L, and c. You could say that I should have let K=IE/L^4/m, and I suppose maybe I should.

4. Thanks to Carl for answering that.

5. Notice that sol1 contains two arbitrary constants _C1 and _C2. In order to plot anything they have to be given concrete values. The two values were taken out of thin air. There is no good reason why you shouldn't choose _C1=1,_C2=0 instead if you are interested in showing the eigen-vibrations (or whatever they are called) corresponding to the eigenfrequencies r[i].
After all a sum like a*cos(w*t)+b*sin(w*t) may be written as A*sin(w*t+phi), and in this context do we care about a phase shift phi or whether the amplitude is A or 1?

Just a comment.
I now remember that I raised the issue of solving for functions of floats in MaplePrimes at

http://www.mapleprimes.com/questions/202411-Solving-For-Function-Of-Float

The problem is illustrated here:

restart;
solve(f(0.1)=5,f(0.1)); #returns NULL
solve(f(1/10)=5,f(1/10)); #OK

To solve the problem we could either convert all floats to rationals by
vals:=convert(vals,rational);

or just convert the arguments to the functions, i.e. the independent variables:

vals:=evalindets(vals,name=float,convert,rational):

thereby keeping the function values as floats.

The reason why my version was slow seems basically to be that in my last loops the procedure D is called Ns*NP*NC times. D is also called when using solve. There is no reason why D should be called in those cases. D should just be replaced with something inert, i.e. just a name like DD. I tried doing it here:
df_dz[s][i,j]:=subs(D=DD,convert(df_dz[s][i,j],D));

Also I skipped
res1:=evalindets(res,function,x->op(0,x)):
and then the last loops would be

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(vals,res,DD[j-NP](z[s][i])(seq(z[s][k],k=NP+1..NP+NC)));
od: od: od:

On a test with (Ns,NP,NC) = (5,50,10) this actually made my version 20% faster than yours in the float case (with the correction for the solve problem mentioned above).

@Carl Love You surely meant: The highest derivative of P in your ODEs is the first. Thus, you can only have three initial conditions, and you can't use D(P) as one of them.

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.