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

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

You will have a better chance of getting a response if you provide the code. You could upload a worksheet.

@samiyare You could try minimizing the sum of squares of pC and pT using LSSolve from the Optimization package.

Input something like

Optimization:-LSSolve([pC,pT]);

It worked for me in Maple 12, which is what I have available this week.

@samiyare You could try minimizing the sum of squares of pC and pT using LSSolve from the Optimization package.

Input something like

Optimization:-LSSolve([pC,pT]);

It worked for me in Maple 12, which is what I have available this week.

@samiyare The algorithm in the link you give could perhaps be adjusted to the present situation like this:

solpar:=dsolve(subs(f(x)=F(x),{eq2=0,eq3=0, C(0)=1,D(C)(0)=c1, T(0)=1, D(T)(0) = t1}), parameters=[c1,t1],numeric,output=listprocedure,known=F);
Cn,Tn:=op(subs(solpar,[C(x),T(x)]));
#Making procedure pC and pT which compute C(b) and T(b), i.e. Cn(b) and Tn(b)
pC:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Cn(b)
end proc;

pT:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Tn(b)
end proc;

#Testing on previously found
sol(0);
pC(-2.11068134505884,-.922230659910390);
pT(-2.11068134505884,-.922230659910390);

Then you would have to use fsolve on the system {pC(c1,t1),pT(c1,t1)}.

@samiyare The algorithm in the link you give could perhaps be adjusted to the present situation like this:

solpar:=dsolve(subs(f(x)=F(x),{eq2=0,eq3=0, C(0)=1,D(C)(0)=c1, T(0)=1, D(T)(0) = t1}), parameters=[c1,t1],numeric,output=listprocedure,known=F);
Cn,Tn:=op(subs(solpar,[C(x),T(x)]));
#Making procedure pC and pT which compute C(b) and T(b), i.e. Cn(b) and Tn(b)
pC:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Cn(b)
end proc;

pT:=proc(c1,t1) if not type([c1,t1],list(numeric)) then return 'procname(_passed)' end if;
solpar(parameters=[c1,t1]);
Tn(b)
end proc;

#Testing on previously found
sol(0);
pC(-2.11068134505884,-.922230659910390);
pT(-2.11068134505884,-.922230659910390);

Then you would have to use fsolve on the system {pC(c1,t1),pT(c1,t1)}.