@pagan Well, I think this is what I meant.
However, I meant not F[i]=sol(y), but F[i]=f(sol(y)).
I hope this will not create any difference.
And yes, without rmemebr option this is rather slow.
I am afraid this will kill the idea, since in reality I need to implement this to a system of 6 DE's...
we will see.
Anyway, this is the closest thing I've got since the beginning of asking questions here.
I guess this should be like
sys:= diff(y(t),t) = cos( 1/2*y(t)-x(t)-Pi^2 ):
for i from 1 to N do
sol[i]:=dsolve([newsys, y(0)=10*sqrt(2)/2], numeric,