@rit I'm not sure I understand what you want to accomplish. Certainly in using matrix multiplication addition is performed, but that would (luckily) not be affected by the device I gave. If it did quite a lot of stuff in Maple would similarly be affected.
A simple example:
restart;
local `+`;
`+`:=proc(a,b) :-`+`(a^~2,b^~2) end proc;
mtaylor(exp(x),x=0,2);
                             1 + x
%;
                             1 + x
eval(%,x=2);
                               3
1+x;
                              2    
                             x  + 1
eval(%,x=2);
                               5

Basically invfourier cannot do the job. invfourier calls fourier and that cannot do it. That it returns 0 is not optimal of course.
I tried replacing all floats by 1:
u:=evalindets(vout_fourier_num2,float,1);
#Then you are dealing with
#u := fourier(phi[3](t), t, w) = w^2*cosh(w)*exp(1+w*(I+w))/(w^4+I*w^3+w^2+I*w+1)
#If you try
invfourier(u,w,t);
# or
fourier(u,w,t);
# you will see that they return partially unevaluated.

@Carl Love I tried some (unsophisticated) code I had written recently. It returned a solution which didn't resemble the one obtained by continuation in lambda. That made me suspicious. Then I recalled reading about the equation of Blasius (and Falkner–Skan) not having unique solutions for certain values of the parameter (Philip Hartman: ODEs, Harold Davis, Intro. to Nonlinear Diff. and Int. Equations, Wikipedia). The present equation is from the same area. 
Hartman imposes the extra condition that the derivative be positive and bounded (by 1) for all t (~eta). A similar condition here would weed out some solutions, but even then it is not clear to me if there is only one. Obviously in Hartman inf = infinity, not 10.

An extra problem here is that there is an apparent singularity at values for theta where 1 - f(eta)^2*lambda = 0, since that quantity appears as the denominator when isolating for f'''.

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