Hi vuishl,

Restart is probably not the best solution, since it only works at the top level (see ?restart) - that is, not inside a loop. But there's a framework that you can use to run these sorts of things much more efficiently, if the objective is to plot: ?dsolve/numeric/Events, which uses the *numeric* ODE solver instead of the *exact* one.

I'm trying to reverse engineer what you're doing with the RootFinding:-NextZero etc., and I think this is more or less what happens: you want to follow a solution until *x*(*t*) is plus or minus one, then negate the derivative of the *x*(*t*). Is that correct? Then a more efficient way to generate such a series of plots starts as follows, in a fresh session:

ODE := diff(diff(x(t), t), t)+.2*(diff(x(t), t))-x(t) = aa*sin(t);

ICs := x(0)=mm, D(x)(0)=nn;

solution := dsolve([ODE, ICs], numeric, parameters=[aa, mm, nn],

events=[[x(t)^2-1, diff(x(t), t)=-diff(x(t), t)]], range=0..100):

Now *solution* is a procedure that contains everything necessary for simulation of your system, including the "bouncing back" of the derivative. The parameters aa, mm, and nn are placeholders that you need to substitute values for before simulation starts (I'll show how in a second). Then you can request the values at a particular point in time, say *t* = 5, as follows:

solution(parameters = [aa = 1.5, mm = 0.999, nn = -2.5]);

solution(5);

This will return a list of equations saying that (approximately) *t* = 5, *x*(*t*) = -0.56, and *diff*(*x*(*t*), *t*) = -1.09. But much more interesting is a routine that will do the phase plot for you:

plots:-odeplot(solution, [x(t), diff(x(t), t)], t=0..100, refine = 1);

This we can use in a loop.

for a from 1.0 to 2.0 by 0.4 do

for m from 0.990 to 1 by 0.004 do

for n from -3.0 to -2.0 by 0.4 do

solution(parameters = [aa = a, mm = m, nn = n]);

plotsetup(ps, plotoutput = sprintf("c:/tmp/partA[%a]M[%a]N[%a].eps", a, m, n), plotoptions = `color,portrait,noborder,width=16cm,height=12cm`);

plots:-odeplot(solution, [x(t), diff(x(t), t)], t=0..100, refine=1);

plotsetup(inline);

end do;

end do;

end do:

As you see, there is no call to *dsolve* inside the loop - that's the beauty of using parameters! All that happens is that the simulation is run and the plot is drawn. Note that the name of the loop variable is different from the name of the parameter - that's always a good idea.

I would try to see if you can run this loop on your older system first. If that doesn't work, we can still break it up into smaller pieces, but I'd think that this is going to be much more efficient than what you had before.

One note: your original version always counted to 200 "bounces", whereas I bounded the simulation by letting it always run to *t*=100 irrespective of the number of bounces. I think we *can* make the simulation always run to 200 bounces, if you'd strongly prefer that; it would be a bit more complicated though. And we'd need a reasonable upper bound on how long that's going to take. Let us know if that is necessary.

Hope this helps,

Erik Postma

Maplesoft.