@Sn4ky Your use of notation for the parameters (EI and L in particular) suggests to me that your pde comes from some branch of physics or engineering. Thus the initial and boundary conditions must be taken from the problem as it appears there.
Were it just a mathematical problem you could play with anything your imagination can come up with.
Not knowing the context, I just selected the conditions using my imagination.

@H-R evalindets may help. Take this (contrived) example:

restart;
S:=Sum(A[n]*cos(n*Pi*x/L)+B[n]*sin(n*Pi*x/L),n=1..N)+Sum(C[n]*cos(n*x/L)+G[n]*sin(n*x/L),n=1..N);
u:=subs(Sin=sin,combine(Sin(m*Pi*x/L)*S)); #To avoid trig simplifications if desired
evalindets(u,specfunc(anything,Sum),s->int(op(1,s),x=-L..L));

@lolly In fact they are different, but not all that much:

resNE:-value()(.345,1);
    [x = 0.345, t = 1., u(x, t) = 9.959477403661344]
resNBE:-value()(.345,1);
    [x = 0.345, t = 1., u(x, t) = 9.989010869640284]

To see the difference in a graph you could do:
p:=(xx,tt)->subs(resNE:-value()(xx,tt),u(x,t))-subs(resNBE:-value()(xx,tt),u(x,t));
#Plotting the difference at time = 1:
plot('p(x,1)',x=0..1,adaptive=false,numpoints=6);

@lolly Now the file upload succeeded.
If you haven't already done so I would suggest that you compare your Neumann version with the one obtained from Maple's pdsolve (the method used is not important for that).
restart;
pde:=diff(u(x,t),t)=(x+t)*diff(u(x,t),x,x)+x+t;
resNE:=pdsolve(pde,{u(x,0)=20*x,D[1](u)(0,t)=0,D[1](u)(1,t)=2*t},numeric,timestep=1/200,spacestep=1/5,method=Euler);
resNE:-animate(t=1);
resNBE:=pdsolve(pde,{u(x,0)=20*x,D[1](u)(0,t)=0,D[1](u)(1,t)=2*t},numeric,timestep=1/200,spacestep=1/5,method=BackwardEuler);
resNBE:-animate(t=1);
#The Dirichlet result you got resembles the Maple result much better.
resD:=pdsolve(pde,{u(x,0)=20*x,u(0,t)=0,u(1,t)=0},numeric,timestep=1/200,spacestep=1/5,method=BackwardEuler);
resD:-animate(t=1);
####################
I took a brief look at your Dirichlet implicit version.
I don't see the point of the double loop:

M is redefined at each iteration and just end up being what corresponds to evaluaing
Matrix(n-1, shape = identity)+a*tau*A/h^2-b*tau
at values for a and b corresponding to the last values of j and m.

Actually, the "exact" result is wrong!
Try
plot(subs(xz_sol,y_sol,[x(t),y(t),z(t)]),t=0..1); p1:=%:
And compare with the corresponding plot from the numerical computation:
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,z(t)]],0..1); p2:=%:
plots:-display(p1,p2);
Indeed it was surprising that the output for y(t) was not a piecewise defined function as it ought to be.

@lolly No, it doesn't. But maybe it is a MaplePrimes problem. To test if there is any problem uploading I shall try here:

test.mw

There was no problem for me.

@lolly Do the links you just posted work when you try yourself? I still get "Bad Request".

None of your links work for me. I get "Bad Request".
In the code posted (as an image, not a good idea) you definitely should not have assignment (:=) but equality (=).

@siamak taghavi You could do like this:

VR22:=0.1178*diff(phi(x),x,x,x,x)-0.2167*diff(phi(x),x,x)+0.0156*diff(psi(x),x,x)+0.2852*phi(x)+0.0804*psi(x);
VS22:=0.3668*diff(psi(x),x,x)-0.0156*diff(phi(x),x,x)-0.8043*psi(x)-0.80400*phi(x);
bok:=evalf(dsolve({VR22=0,VS22=0}));

PHI,PSI:=op(subs(bok,[phi(x),psi(x)]));
Eqs:={eval(PHI,x=.83)=1,eval(diff(PHI,x),x=0.83)=0,eval(PHI,x=-.83)=1,eval(diff(PHI,x),x=-0.83)=0,
eval(PSI,x=.83)=1,eval(PSI,x=-0.83)=1};
C:=fsolve(Eqs,indets(%,name));
eval(bok,C);
SOL:=fnormal(evalc(%));
plot(rhs~(SOL),x=-.83..0.83);

@Markiyan Hirnyk There are (at least) 2 solutions for Q. For N=22 fsolve(Q,0) gave the negative one.
With fsolve(Q,0.002) we get the positive one.
I added a remark to my answer above.

Runge-Kutta methods are used for solving initial value problems. You have a boundary value problem.
So please elaborate.

@Preben Alsholm I tried the continuation I mentioned above, but had only success with A1=4.5.
So just a slight improvement from A=0.43 obtained without continuation.
Then I tried using the result obtained for A=0.45 by continuation from A=0.2 as an aproximate solution for the problem with A=0.5. It didn't work. But I was unpleasantly surprised to find that the continuation solution didn't even work when A=0.45! That is: I'm handing dsolve the solution, but yet it cannot find the solution with that as a start!!

sys:=eval({Eq1, Eq2, Eq3, Eq4, Eq5, Eq6, bcs1},A=0.2*(1-c)+c*0.45);
res:=dsolve(sys, [f(eta), F(eta), G(eta), H(eta), theta(eta), theta1(eta)], numeric, output = listprocedure,continuation=c); #Works
plots:-odeplot(res,[[eta,diff(f(eta),eta)],[eta,F(eta)]],0..blt,linestyle=3);
#Now dsolve is given sys with A=0.45 the approximate solution (res) found above for A = 0.45 by continuation from A=0.2, yet we get the usual error:
res2:=dsolve(eval({Eq1, Eq2, Eq3, Eq4,Eq5,Eq6,bcs1}, A = 0.45),numeric, output = listprocedure,approxsoln=res); #Does not work!!

 

spline is deprecated, use Spline from the CurveFitting package.
Maybe you could use CurveFitting:-ArrayInterpolation?

@cloz54 Cutting and pasting your module definition into my version of Notepad and saving it as a textfile. Reading it into Maple 18 I got

                     "Est?n en perspectiva"
                      "Spain means Espa?"
Actually the question marks appeared as squares.

I also tried executing the moduledefinition in Maple, saving it using
save TestModule,"F:/MapleDiverse/MaplePrimes/testmod.mpl";

and then reading it into a fresh worksheet:

read("F:/MapleDiverse/MaplePrimes/testmod.mpl");

Then I got:
with(TestModule):
test();
                     "Están en perspectiva"
                      "Spain means España"
                               0

@Axel Vogt I tried the following. All had the same result on my computer:

print("Están en perspectiva");
p:=proc() print("Están en perspectiva") end proc;
p();
m:=module() export p; p:=proc() print("Están en perspectiva") end proc; end module;
m:-p();
Q:=proc() module() export p; p:=proc() print("Están en perspectiva") end proc; end module end proc;
m:=Q();
m:-p();