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

It is not clear to me what you are trying to do. What role does the variable 't' play?
Maybe you could tell us where this problem came up?

@samiyare Notice that I have added a faster version in my answer above.