@JJames _Z1 (printed with a trailing tilde meaning assumptions made on _Z1) stands for an arbitrary integer. So if 'res' above contains _Z1 you can replace it with any integer e.g. 1. 'res' actually also contains a name _C3 (without a trailing tilde). A name like that stands for an arbitrary complex constant.

eval(res,{_Z1=1,_C3=I/2});

As in solving the eigenvalue problem for a matrix A: Find the numbers lambda s.t. Av = lambda*v has a nontrivial solution, there may be ways of finding lambda before actually finding the corresponding v's. In this case we solve the characteristic equation det(A-lambda*Id) = 0.

Above we found the possible values for lambda (actually assuming that they were positive). There were infinitely many, and I chose to use a simpler name n instead of _Z1 or _Z2 or what have you (Maple changes the number if you do it again without a restart).

lambda = (a^2+Pi^2*n^2)^(3/2)/a.

After having found the possible values of lambda you can find all the corresponding eigenfunctions:

dsolve(eval({ode1,ode2,bcs},lambda=(a^2+n^2*Pi^2)^(3/2)/a)) assuming n::integer;

However, Maple can (in principle) solve the whole problem in one command as I noticed in the corrected version of my first answer to your question:

dsolve({ode1,ode2,bcs},{alpha(z),H(z),lambda});

The output, however, is rather unwieldy, but hopefully correct.

With a little more work it can be seen that the corresponding sequence of negative values
lambda= - (a^2+n^2*Pi^2)^(3/2)/a)
are also eigenvalues (n<>0 though!).

@JohnJames There is no ready made procedure, so yes, you will hve to do it manually. The individual manual steps could, however, be done in Maple, but you will most likely find it much easier to do them on paper.

To illustrate that last point I can give you a few steps, which will probably convince you:

res1:=dsolve(odes[1]) assuming _c[1]<0; #if that is the relevant assumption
res2:=dsolve(odes[2]);
eval(IBC,f=unapply(rhs(op(1,res)),x1,x2));
#We want a nontrivial solution, so
eval(%,{_F1(x1)=1,_F2(x2)=1});
diff(res1,x1);
eval(%,x1=0);
eval(rhs(%%),{x1=L1,_C1=0});
solve(%=0,{_c[1]},AllSolutions);

etc.

In other words, you will have to know what to do yourself at each step. Doing it in Maple just adds to the difficulty in this case.

@JohnJames There is no ready made procedure, so yes, you will hve to do it manually. The individual manual steps could, however, be done in Maple, but you will most likely find it much easier to do them on paper.

To illustrate that last point I can give you a few steps, which will probably convince you:

res1:=dsolve(odes[1]) assuming _c[1]<0; #if that is the relevant assumption
res2:=dsolve(odes[2]);
eval(IBC,f=unapply(rhs(op(1,res)),x1,x2));
#We want a nontrivial solution, so
eval(%,{_F1(x1)=1,_F2(x2)=1});
diff(res1,x1);
eval(%,x1=0);
eval(rhs(%%),{x1=L1,_C1=0});
solve(%=0,{_c[1]},AllSolutions);

etc.

In other words, you will have to know what to do yourself at each step. Doing it in Maple just adds to the difficulty in this case.

Relying on the order and form of the output from convertsys you could define vn by using op:

vn:=op([2,1,2,1],Deq);

As for eval(%,dat[2]); I don't know what your intention is (in view of the next line with remove).

Relying on the order and form of the output from convertsys you could define vn by using op:

vn:=op([2,1,2,1],Deq);

As for eval(%,dat[2]); I don't know what your intention is (in view of the next line with remove).

@herclau Yes, I have experienced several times that copying from a worksheet to MaplePrimes fails to include everything. This, unfortunately, makes it necessary to read the whole code after copying, something that I clearly failed to do with Dsolve2.

@herclau Yes, I have experienced several times that copying from a worksheet to MaplePrimes fails to include everything. This, unfortunately, makes it necessary to read the whole code after copying, something that I clearly failed to do with Dsolve2.

@herclau
I don't get any email notifications either anymore!

Whereas LinearAlgebra is a module, DEtools is not, which explains the difference in their treatment below.

I have made a few changes in your procedure besides making it take different input (and output).

Dsolve2:=proc (eqs::{list,set}(equation),ics::{list,set}(equation),vars::{set,list}(function(name)))
local t,dat,A,b,v,vi,ic,Xt,X,x,sys;
uses LinearAlgebra;
t:=op(op~(vars));
dat:=DEtools[convertsys](eqs,ics,vars,t,x);
A,b:=GenerateMatrix(rhs~(dat[1]),lhs~(dat[2]));
Xt:=applyop~(apply,1,dat[2],t);
X:=lhs~();
v:=diff~(X,t)=A.X-b;
vi:=eval(X,t=dat[3])=;
sys:=convert(Equate(lhs(v),rhs(v)),set) union convert(Equate(lhs(vi),rhs(vi)),set);
userinfo(2,'procname',print(v,vi,));
eval(dsolve(sys,_rest),Xt);
end proc:

infolevel[Dsolve2]:=2:
Dsolve2({deq},{u(0)=2,D(u)(0)=-1},{u(t)});
Dsolve2({deq},{u(0)=2,D(u)(0)=-1},{u(t)},numeric,output=Array([0,0.25,0.5,0.75,1]));

#A comment: Your procedure as well as my modified version obviously only handles linear systems.

@herclau
I don't get any email notifications either anymore!

Whereas LinearAlgebra is a module, DEtools is not, which explains the difference in their treatment below.

I have made a few changes in your procedure besides making it take different input (and output).

Dsolve2:=proc (eqs::{list,set}(equation),ics::{list,set}(equation),vars::{set,list}(function(name)))
local t,dat,A,b,v,vi,ic,Xt,X,x,sys;
uses LinearAlgebra;
t:=op(op~(vars));
dat:=DEtools[convertsys](eqs,ics,vars,t,x);
A,b:=GenerateMatrix(rhs~(dat[1]),lhs~(dat[2]));
Xt:=applyop~(apply,1,dat[2],t);
X:=lhs~();
v:=diff~(X,t)=A.X-b;
vi:=eval(X,t=dat[3])=;
sys:=convert(Equate(lhs(v),rhs(v)),set) union convert(Equate(lhs(vi),rhs(vi)),set);
userinfo(2,'procname',print(v,vi,));
eval(dsolve(sys,_rest),Xt);
end proc:

infolevel[Dsolve2]:=2:
Dsolve2({deq},{u(0)=2,D(u)(0)=-1},{u(t)});
Dsolve2({deq},{u(0)=2,D(u)(0)=-1},{u(t)},numeric,output=Array([0,0.25,0.5,0.75,1]));

#A comment: Your procedure as well as my modified version obviously only handles linear systems.

There is no infinite series, so there is no question of convergence:
The idea in Gantmacher is to replace f by a polynomial r which has the "same values on the spectrum" as f. By the same values is meant not only r(L) = f(L) for every eigenvalue but also equality for the derivatives up to one less than their individual multiplicities in the minimal polynomial of the matrix.
With your new matrix:
G := Matrix(3, 3, {(1, 1) = 218/29, (1, 2) = 96/29, (1, 3) = 0, (2, 1) = 120/29, (2, 2) = 130/29, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 5})

Eigenvalues(G);
MinimalPolynomial(G,lambda); #in this case (as for J) the same as the characteristic poly
factor(%);
#The 3 multiplicities for 2, 5, and 10 are 1:
#m[L]=1, for L=2,5,10. Thus m[L]-1 = 0 for all three (see eqs below).
r:=unapply(add(a[k]*lambda^k,k=0..2),lambda);
f:=x->ln(1+x);
#The set of equations in this case. In the inner seq k runs from 0 to m[L]-1 (in this case 0):
eqs:=op~({seq(seq([(D@@k)(r)(L)=(D@@k)(f)(L)],k=0..),L=[2,5,10])});
solve(%);
R:=eval(r(lambda),%);
M:=eval(R,lambda=G);
#Check:
MatrixExponential(M)-1;

There is no infinite series, so there is no question of convergence:
The idea in Gantmacher is to replace f by a polynomial r which has the "same values on the spectrum" as f. By the same values is meant not only r(L) = f(L) for every eigenvalue but also equality for the derivatives up to one less than their individual multiplicities in the minimal polynomial of the matrix.
With your new matrix:
G := Matrix(3, 3, {(1, 1) = 218/29, (1, 2) = 96/29, (1, 3) = 0, (2, 1) = 120/29, (2, 2) = 130/29, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 5})

Eigenvalues(G);
MinimalPolynomial(G,lambda); #in this case (as for J) the same as the characteristic poly
factor(%);
#The 3 multiplicities for 2, 5, and 10 are 1:
#m[L]=1, for L=2,5,10. Thus m[L]-1 = 0 for all three (see eqs below).
r:=unapply(add(a[k]*lambda^k,k=0..2),lambda);
f:=x->ln(1+x);
#The set of equations in this case. In the inner seq k runs from 0 to m[L]-1 (in this case 0):
eqs:=op~({seq(seq([(D@@k)(r)(L)=(D@@k)(f)(L)],k=0..),L=[2,5,10])});
solve(%);
R:=eval(r(lambda),%);
M:=eval(R,lambda=G);
#Check:
MatrixExponential(M)-1;

@Markiyan Hirnyk I didn't mean to imply equivalence. But from the differential inequation follows the integrated version.

@Markiyan Hirnyk I didn't mean to imply equivalence. But from the differential inequation follows the integrated version.

@Markiyan Hirnyk If N(t) is your proposed

N1:=t->-(-Pi+Pi*exp(-t))/exp(-t)-2-sin(exp(t^2));

then N1(0) = -2-sin(1). But you assumed N(0) = 0.

@Markiyan Hirnyk If N(t) is your proposed

N1:=t->-(-Pi+Pi*exp(-t))/exp(-t)-2-sin(exp(t^2));

then N1(0) = -2-sin(1). But you assumed N(0) = 0.