The answer to this is very similar to the answer to your last question. If you use assumptions and simplify the dsolve results with evalc, then those results will simplify dramatically. Then you can check directly that they satisfy the ODEs. (I am not saying that there's anything wrong with odetest or Mehdi's Answer; I am showing that Maple can work well with imaginary values.)

restart:
assume(a > 0, -Pi < m, m <= Pi);
sys_ode:=
     diff(d11(m),m) =
           -(3*sin(m)^2-1)*d31(m)/a^(3/2)+(-3*cos(m)*sin(m)/a^(3/2))*d41(m),
     diff(d21(m),m) =
          (-3*cos(m)*sin(m)/a^(3/2))*d31(m)-(3*cos(m)^2-1)*d41(m)/a^(3/2),
     diff(d31(m),m) = -a^(3/2)*d11(m),
     diff(d41(m), m) = -a^(3/2)*d21(m)
:
Sol1:= dsolve({sys_ode}):
Sol2:= simplify(evalc(Sol1));

Now verify that these solutions satisfy the original system of ODEs:

simplify((lhs-rhs)~(eval({sys_ode}, Sol2)));

There are very many ways to do this in Maple. Here's a way that puts the values in a Vector:

f:= x-> x^3+2*x+1:
fvalues:= Vector(10001, k-> f((k-1)*0.001));

What do you want to do with the values? Knowing that makes it easier to decide which way to generate them.

You forgot the parentheses. Surrounding a block with (*  *) comments it out. Also, # comments until the end of the line.

Use numeric integration: Enter Int with a capital I, and then

evalf(%);

0.1780721092

Such an operation is not usually done with a for loop in Maple, although it can be done that way. It is done like this:

a:= proc(n::nonnegint)
option remember;
local t;
     unapply(f(thisproc(n-1)(t), thisproc(n-2)(t)), t)
end proc:
a(0):= 0:  a(1):= t-> t:

When you say that you want to plot a[N](t), do you mean for a specific N? Note that with the procedure a as defined above, the notation is a(N)(t) rather than a[N](t).

The first error is that you define the procedure named GenerateStencil, but in PoissonSolve you call Stencil, not GenerateStencil.

There are other errors, but that's a start.

It is -m, but modulo 2*Pi, depending on the branch of the logarithm.

simplify(evalc(I*ln(cos(m)+I*sin(m)))) assuming -Pi < m, m <= Pi;

                                                           -m

Is the assumption -Pi < m <= Pi valid in your case?

First, end each of the EQs with a semicolon. Then assign the initial conditions to a variable, say ics, and end with a semicolon. Then

dsolve({EQ||(1..5), ics}, {q||(1..5)(T)}, method= laplace, convert_to_exact= false):
allvalues(%):  #obtain approximate roots of polynomials.
evalc(%):  #convert exp(a*I*T) to sin cos forms.
simplify(%, zero);  #Remove spurious imaginary parts.

a:= proc(n) option remember; a(0)+a(n-1) end proc:  a(0):= 1:
a(4);

                                                  5

Your phi46N in the matrix is what is known as an "escaped local" variable. (This is the problem that you suspected.) So, it is different from the variable phi46N which you enter at the top level. To refer to this variable, you need to extract it from the context in which it occurs. So one form of a correct solve command is

solve({seq(eq[i], i= 1..6)}, convert(phi4N, set));

Constructing the eq[i] in a loop is actually superfluous, so you could do

solve(convert(phiN =~ phi4N, set), convert(phi4N, set));

Set Digits to a higher value (and higher than 15), redo the computation, and subtract the new results for the eigenvalues from the original. You can do the same with the eigenvectors except that they made be reversed in sign.

When the number variables being solved for is less than the number of equations you should generally use eliminate instead of solve.

eliminate(sys, drivers);

Note that eliminate always returns pairs of sets: The first set is the elimination equations; the second set is the original equations with the appropriate variables eliminated via the substitutions from the first set.

In the pdsolve command, use the options timestep and spacestep. For example,

pds := pdsolve(pde, IBC, numeric, timestep=.01, spacestep= 1);

I cannot vouch for the accuracy of the plots that you get with this, but you will get the plots instead of error messages. I tried smaller values, but the left sides of the plots still have that sharp drop. But note that the boundary condition u(t,-50) = 0 necessitates that the left sides go to 0, so maybe those plots are accurate.

Here's how to make the input the differential equation itself. The procedure preprocessODE takes the ODE (given in the standard Maple form accepted by dsolve) and returns the function f that is passed to eul. I made a slight modification to eul so that it calls preprocessODE. (I am not sure if this will work in Maple 5; the error command syntax might be different. Let me know.)

preprocessODE:= proc(eq::`=`)
local d, y, x, f;
    d:= indets(eq, specfunc(anything, diff));
    if nops(d) = 0 then error "No derivative found." end if;
    if nops(d) > 1 then error "More than 1 derivative found." end if;
    d:= d[];
    y:= op([1,0], d);  x:= op([1,1], d);
    f:= {solve(eq, d)};
    if nops(f) <> 1 then error "Cannot isolate derivative." end if;
    unapply(subs(y(x)= y, f[]), [x,y])
end proc:
    
eul:=proc(ODE::`=`,h,x0,y0,xn)
  local f,no_points,x1,y1,i;
  f:= preprocessODE(ODE);
  no_points:=round(evalf((xn-x0)/h));
  x1:=x0;
  y1:=y0;
 
  for i from 1 to no_points do
      y1:=y1+evalf(h*f(x1,y1));
      x1:=x1+h
  end do;
  y1
end proc:
       

Example of use:

eul(diff(y(x), x) = x^2*y(x)^3, .01, 0, 1, 1);

1.68181553277060

You cannot use indexed variables where the index or the base are the same as the variable of integration. So you need to change lambda[v], lambda[t], and t[c]. I changed them to L__v, L__t, and t__c  (which print as subscripted variables in Maple 17 & 18). Also, you need to assume that the lambdas are positive.

int(
     int(
          L__v*L__t*exp(-L__v*v-L__t*t),
          v = (1/2)*(q[p]+q[p]*t__c*t+2*S[di]*h*t)/(h*t) .. infinity
     ), t = 0 .. infinity
) assuming L__v > 0, L__t > 0;