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.