@Carl Love I just tried giving a vote up. Now the OP has a reputation of 5.

@marc sancandi I admit to have been put off by your heavy use of sarcasm, which I think is not conducive to a sober and rational discussion.

Does the problem you have show up in an example simple enough to view on a computer screen?
If so, why not upload a worksheet not relying om an .m file?

In your reply to Tom Leslie you say that the system is huge, so I guess that for now you don't have a simpler example.
Often enough I have heard (or read) numerical analysts say that for any numerical method there exist inputs for which the method won't work.

Should we just trust your statement:

# File "MC.m" below contains an example of one of them (a second degree ODE
# reduced to a couple of 2 first order ODEs ; details about these equations
# are of no importane here).

?

@shahid Try not to assign to RE and then do:

res := dsolve(eval({bc, eq1, eq2},RE=15), numeric); #Works for me

Now you can experiment with continuation in RE or with using an approximate solution to get success for higher values of RE.
See  ?dsolve,numeric,bvp

@xhimi I make 4 (randomly chosen) matrices with size 16x16 and with complex entries (no symbolic entries). I give them shape=hermitian; then I know that the eigenvalues will be real and more importantly, Maple knows it too.
If you happen to have an hermitian matrix and the matrix is not given the shape hermitian, then there will no doubt be imaginary roundoff errors in the eigenvalues. Those can be removed by the use of fnormal and then simplify. For an example see the bottom.
restart;
with(LinearAlgebra):
n:=16:
p:=RandomTools:-Generate(complex(float),makeproc);
A:=RandomMatrix(n,generator=p,shape=hermitian);
B:=RandomMatrix(n,generator=p,shape=hermitian);
C:=RandomMatrix(n,generator=p,shape=hermitian);
D1:=RandomMatrix(n,generator=p,shape=hermitian);
L:=[A,B,C,D1];
LE:=[Eigenvectors]~(L):
LEV:=map2(op,2,LE);
LEV[1]; #Matrix of eigenvectors of A
seq(LEV[1][..,i],i=1..n); #Eigenvectors for A
m:=nops(L);
res:=IntersectionBasis([seq([seq(LEV[j][..,i],i=1..n)],j=1..m)]);
interface(rtablesize=16);
res;
LE[1][1]; #Eigenvalues for matrix A
##############################
D2:=Matrix(D1): #Making a new matrix with the same elements as D1, but with default options
MatrixOptions(D2);
ev:=Eigenvectors(D2):
ev[1]; #Vector of eigenvalues
simplify(fnormal(ev[1])); #Stripping roundoff errors

@shahid You have a new equation eq1.
Secondly, you cannot use [] instead of parentheses ().

Change that before attempting again. You may have to work harder after that.

@xhimi To see matrices or vectors with dimensions higher than 10 do:
interface( rtablesize= 16); #
to set it at 16. By using
interface( rtablesize= infinity);
all sizes will be shown.
###
You say that your code stops at m := nops(L);
The answer obviously should be 4 in your case.
##
What type of elements do your matrices have? If they have datatype=float there ought not be any problem.
If on the other hand your entries have symbolic values like names or exact expressions like sqrt(2) or sin(1) then you can easily run into problems.
Remember that to find the eigenvalues (and the eigenvectors) you have to find the roots of a polynomial of degree 16, which is impossible in general except numerically.

@Carl Love Well, you take the real part. The real part does become zero:

evalc(ex1(x,t,z));
          -1.132*10^11*exp(9.9*10^6*x)*cos(-1.95*10^6*z+2.98*10^15*t)

Thus it is zero where cos(-1.95*10^6*z+2.98*10^15*t)=0, i.e. for
-1.95000000*10^6*z+2.980000000*10^15*t = Pi/2+p*Pi, p any integer and for any x.

## Getting rid of large numbers makes this less confusing:
ex1:= (x,t,z)-> Re(-exp(x + I*(z-t)));
plots:-implicitplot3d(
     ex1(x,t,z), x= -10..0, t= 0..10, z= 0..10,
     axes= boxed, style= patchcontour, scaling= constrained, shading= z,
     grid= [20$3]
);

@xhimi Well, here is a way. I still make the matrices small so that results are easy to understand.

restart;
with(LinearAlgebra):
n:=5: #Using nxn matrices
## Taking just 3 matrices:
A:=RandomMatrix(n,datatype=float);
B:=RandomMatrix(n,datatype=float);
C:=RandomMatrix(n,datatype=float);
L:=[A,B,C]; #Putting the matrices in a list
LE:=[Eigenvectors]~(L): #Finding vectors of eigenvalues and corresponding matrices whose columns are eigenvectors
LE[1]; #Result for A
LEV:=map2(op,2,LE); #Leaving out the vectors of eigenvalues
LEV[1]; #Matrix of eigenvectors of A
seq(LEV[1][..,i],i=1..n); #Eigenvectors for A
m:=nops(L); #Number of matrices
IntersectionBasis([seq([seq(LEV[j][..,i],i=1..n)],j=1..m)]);

xhimi 0

You say that you have a bunch of matrices. Why not just find all the eigenspaces and then use IntersectionBasis in one call?

@Preben Alsholm  I tried to move your reply to my answer to where it belongs, but it appears as my reply.

xhimi 0  WROTE:

Thanks a lot!

 I have tried the intersection basis but I am not not able to have a procedure that does both processes at the same time, i.e., finding the eigenspace and then recursively intersecting. HELP PLEASE!

Could you give us an actual example? I don't quite understand what you are trying to say.

Obviously the equation LambertW(ln(s)) = ln(z) doesn't have s=1 as a solution generically, i.e unless z=1, since
eval(LambertW(ln(s)) = ln(z), s=1) returns 0=ln(z). So I see that as a bug.
It is strange that occasionally that error occurs.
I tried
for sym in {a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,x,y,z} do solve({LambertW(ln(zz)) = ln(sym)}, zz) end do;
and didn't see any problems, though.
##Added: I also tried:
for sym in {a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,t,u,v,x,y,z} do solve({LambertW(ln(s)) = ln(sym)}, s) end do;
which exhibited the error for sym= t,u,v,x,y,z.
It so happens that t,u,v,x,y,z all follow after s in the alphabet.
## Suspicion confirmed by now lettin g be the unknown:
for sym in {a,b,c,d,e,f,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,x,y,z} do solve({LambertW(ln(g)) = ln(sym)}, g) end do;

Clearly a bug and also prsent in Maple 2016.
I shall submit an SCR.


 

@GoitaHass You are not using the code I gave you in your attached worksheet. Your code contains several syntactical errors. I wasn't blunt enough when I wrote: Assuming that your 2-D input can be interpreted to mean the same as this corrected 1-D version ..

Try to take a copy of my reply above and paste it into a simple text editor like NotePad in Windows. You can make Maple work like that.
Since you are new to Maple you may not be aware of the important difference between 2D input and 1D input (aka Maple input).

Go to the menu Tools/Options/Display/Input display choose Maple Notation.
Keep the Options window open.
After that in the Options window go to the tab Interface/Default format for new worksheets choose Worksheet.
Finally click the button Apply Globally.

This action can always be reversed.

Let me add that my personal choice would be to choose simpler names as in this modified version:
restart;
ode := mu1*E1(x)-mu2*(u(x)-d1*E1(x))/d2+diff(gamma1*E1(x)-gamma2*(u(x)-d1*E1(x))/d2, x) = 0;
isc:=E1(0) = 0;
dsolve({ode,isc},E1(x));

@John Fredsted Will this one work in your situation?

conv:=proc(u) local g, res;
  res:=evalindets(u,`*`,s->g(selectremove(type,s,{identical(tau),identical(f(t))})));
  res:=subs(tau*f(t)=f(-t),res);
  eval(%,g=`*`)
end proc;

conv(expr);
conv(tau*f(t)+6*f(t)-(1+7*I)*tau*f(t));
conv([expr,tau*f(t)+6*f(t)-(1+7*I)*tau*f(t)]);

#But maybe just evalindets with algsubs?

evalindets([expr,tau*f(t)+6*f(t)-(1+7*I)*tau*f(t)],`*`,s->algsubs(tau*f(t)=f(-t),s));

## Well, maybe I'm not quite getting it: You probably have many different functions f?
## That could be handled as in this example:
conv:=proc(u) local g, res, tps,S1,S; global tau,t,f,F,h;
  S1:={f(t),F(t),h(t)}; #Just an example
  tps:=identical~({tau} union S1);
  S1:=convert(S1,list);
  S:=tau*~S1=~eval(S1,t=-t);
  res:=evalindets(u,`*`,s->g(selectremove(type,s,tps)));
  res:=eval(res,S);
  eval(%,g=`*`)
end proc;
### And a version with algsubs combined with foldr and evalindets:
L:=[f(t),F(t),h(t)];
S:=tau*~L=~eval(L,t=-t);
conv:=u->foldr(algsubs,u,op(S));
evalindets([expr,tau*f(t)+6*h(t)-(1+7*I)*tau*F(t)],`*`,conv);