Question: code efficiency

Dear all,

Some time ago I asked a question concerning the algorythm for numerical dsolve and related issues here.

I got a nice answer from pagan and implemented the algorythm.

However, since my system is big, the iterations take long time.

Here is the code itself:

sys1:=[TRIC[1],TRIC[2],TRIC[3],TRIC[7],TRIC[8],TRIC[9]]:
> sys2:=[TRIC[4],TRIC[5],TRIC[6],TRIC[10],TRIC[11],TRIC[12]]:
> GRAD:=subs(u=FF[1](t),a=FF[2](t),k=k(t),tau=tau(t),psi[1]=psi[1](t),psi[2]=psi[2](t),psi[3]=psi[3](t),[diff(HAMN,u),diff(HAMN,a)]);
> REC:=[FF[1](t)+alpha[1]*GRAD[1],FF[2](t)+alpha[2]*GRAD[2]];
> F[0]:=t->0.5:
> G[0]:=t->0.2:
> GR[0]:=0:
> DR[0]:=0:
> KK[0]:=0:
> N:=10:
> alpha[1]:=0.01:
> alpha[2]:=0.01:

for i from 1 to N do
> newsys1:=subs(u(t)=F[i-1](t),a(t)=G[i-1](t),sys1):
> print(newsys1);
>  sol1[i]:=dsolve(newsys1, numeric,
>                   output=listprocedure, known=[F[i-1],G[i-1]]):
> KK[i]:=subsop(3=remember,eval(k(t),sol1[i])):
>
> MM[i]:=subsop(3=remember,eval(m(t),sol1[i])):
> TT[i]:=subsop(3=remember,eval(tau(t),sol1[i])):
>
>
> newsys2:=subs(u=F[i-1],a=G[i-1],k=KK[i],m=MM[i],tau=TT[i],sys2):
> sol2[i]:=dsolve(newsys2, numeric,
>                   output=listprocedure, known=[F[i-1],G[i-1],KK[i],MM[i],TT[i]]):
> F[i]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i]))(t),k(t)=subsop(3=remember,eval(k(t),sol1[i]))(t),
> tau(t)=subsop(3=remember,eval(tau(t),sol1[i]))(t),                              FF[1]=F[i-1],
>                               REC[1])),
>                  t,numeric,proc_options=[remember]):
> print(newsys2);
> G[i]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i]))(t),psi[3](t)=subsop(3=remember,eval(psi[3](t),sol2[i]))(t),k(t)=subsop(3=remember,eval(k(t),sol1[i]))(t),
> tau(t)=subsop(3=remember,eval(tau(t),sol1[i]))(t),                              FF[2]=G[i-1],
>                               REC[2])),
>                  t,numeric,proc_options=[remember]):
> GR[i]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i]))(t),psi[3](t)=subsop(3=remember,eval(psi[3](t),sol2[i]))(t),k(t)=subsop(3=remember,eval(k(t),sol1[i]))(t),tau(t)=subsop(3=remember,eval(tau(t),sol1[i]))(t),FF[1]=F[i-1],FF[2]=G[i-1],
>                               GRAD[1])),
>                  t,numeric,proc_options=[remember]):
>
> DR[i]:=unapply('evalf'(subs(psi[1](t)=subsop(3=remember,eval(psi[1](t),sol2[i]))(t),psi[3](t)=subsop(3=remember,eval(psi[3](t),sol2[i]))(t),k(t)=subsop(3=remember,eval(k(t),sol1[i]))(t),tau(t)=subsop(3=remember,eval(tau(t),sol1[i]))(t),FF[1]=F[i-1],FF[2]=G[i-1],
>                               GRAD[2])),
>                  t,numeric,proc_options=[remember]):
>
> end do:

sys1 and sys2 are two ODE systems with 3 equations. They have to be solved sequentially...

GRAD is a 2-dim vector function.

10 iterations take app. 40 minutes, while I need more than 10.

Are there any ideas as to how I may optimize this to take less time?

Thanks for any help.