@acer Just for the record:

shift:=(f::procedure)->(x->f(x+1)):
> shift(sin);

                            x -> f(x + 1)

> "(Pi/6-1);

                              f(1/6 Pi)

> kernelopts(version);

  Maple V, Release 4, IBM INTEL NT, Dec 15, 1995, WIN-5403-951412-1

#Notice that before the introduction of strings in Maple the ditto was " instead of %.

I have not tried your code (it is not easy to copy and paste with all that output included), but since you are aware of the option maxfun you could simply increase that or set it at 0 (which actually means infinity).

I tried your definition in Maple 18 and Maple 11. I obtained what you wanted. So how come you don't?

@xmmood The reason why dsolve/numeric complains when being fed the DAE system must be that the code for detecting imaginary numbers in the input is rather simple. The problem is that the imaginary number I is seen as appearing as a part of the index in S[b,1,I] and in X[b,1,I] (!!!).
That problem is easily solved.
Starting from your original system, which I called sys, you can do:

eps:=0.001:
sysA:=subs(X[b, 1, S](t)/X[b, 1, B, H](t)=piecewise(X[b,1,B,H](t)<eps,1,X[b, 1, S](t)/X[b, 1, B, H](t)),sys);
res:=dsolve(subs(I=ii,sysA),numeric,maxfun=10^5);
#To verify that "I" is the problem try removing subs (or changing ii to I).
depvar:=subs(I=ii,convert(indets(rhs~(sysA),function(identical(t))),list));
plots:-odeplot(res,[seq([t,log10(depvar[i])],i=1..13)],0..100,legend=[seq(depvar[i],i=1..13)]);
#Instead of ii you can use whatever you like (to some degree) e.g. the name ` I`. Notice the space before I.
###############
I shall report the following simple version as an SCR:
restart;
sys:={diff(x[I](t),t)=2,x[I](0)=0};
res:=dsolve(sys,numeric);
res(3);
sys2:={diff(x[I](t),t)=a(t),x[I](0)=0,a(t)=2};
res2:=dsolve(sys2,numeric);
res2a:=dsolve(subs(I=ii,sys2),numeric);
res2a(3);

Could you show us your system? Please use text so we can copy and paste. Alternatively, you can upload a worksheet.

What is your input into Maple (text please!).

@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);