And so with this provocative title, "pushing dsolve to its limits" I want to share some difficulties I've been having in doing just that. I'm looking at a dynamic system of 3 ODEs. The system has a continuum of stationary points along a line. For each point on the line, there exist a stable (center) manifold, also a line, such that the point may be approached from both directions. However, simulating the converging trajectory has proven difficult.

I have simulated as follows: reverse time and treat the system as an initial value problem. Thus I take initial values for the 3 variables of the system near the stationary point (but not exactly on it) and run time backwards: the system diverges away from the stationary point on the stable manifold, thus describing the trajectory of interest.

The problem is that the system is very stiff. For some points on this continuum (near the extremities) I've been unable to get a satisfactory numerical simulation of the stable trajectory.

This is what the simulated dynamics look like. This is (if it displays properly) an animated gif of 50 trajectories associated with 50 stationary points on the continuum.

[ I will attempt to insert the animated gif in a comment below, as it's not loading properly at this time]

A still picture of these trajectories from one side ("above"):

A still picture of these trajectories from the other side ("below"):

It is clear that the dynamics from above are very nonlinear. As a result, I'm unable to get some of the trajectories to converge as (I think) they should.

A picture with fewer trajectories, which may be easier to read:

Along the trajectory, one of the variables, q(t), changes its behaviour dramatically. In a (very tiny) neighborhood of the stationary point, q(t) is declining at a rate that is detectable by the system, namely from a value that is one epsilon short of 0.03 to a lower value such as 0.029 (depending on which point we are studying) in 0.1 units of time. Away from the stationary point, this same variable q(t) is rising, from values like 0.029 to 0.03 in about 80 units of time.

This behaviour can be seen in the figures below. Thus q(t) changes slope dramatically from very flat to much steeper. For some reason this seems to be difficult for dsolve to detect for some of the points.

I have tried every trick I know. I'm looking for new ideas.

I have set high values of digits, low values of relerr, abserr, initstep, maxstep, I have used the option 'compile", I have experimented with the different methods, to no avail -- what convergence I do get, I get with rk45, with digits set to anything between 8 and 20, and with or without all the fancy error-control options. All the options do is to alter the time it takes to simulate, but the results are always basically the same: convergence from below is reasonably easy to obtain, convergence from above fails for some points on one side of the continuum, i.e. the red-coloured trajectories in the picture.

Any suggestions welcome!