@pagan Hi,

Your suggestion turned to be very useful.

However now I am getting stuck with implementation of the same thing with two different functions.

It seems, that the algorithm although updating both ogf them with every i, do not optimize w.r.t. the second of them , a.

In the following code I use the sequential solution of two systems of DEs (called newsys1 and newsys2) with two functions, F[i], G[i] which are updated according tot he solution received at i-th step. However, the results do not differ much for G[i]=const: the same form of F[i] function i obtained. Is there a mostake in my code, or I have to use some sort of a double.cycle algorithm instead of one here?

Here is the code:

for i from 1 to N do

>

> print(i);

> 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:

which is the version of what You suggested.

Thanks for any hints