###
The Steinmetz-Larter example

In the peroxidase-oxidase reaction a peroxidase catalyzes an aerobic oxidation:

where () is a general electron donor. Under the proper conditions this reaction yields oscillations in the concentrations of , and various reaction intermediates. Steinmetz and Larter [Steinmetz and Larter 1991] derived the following system of four coupled nonlinear differential rate equations:

| (32)

where are state variables and , , , , , , , , and are parameters. Although highly simplified, the model is able to reproduce the three modes of simple, chaotic, and bursting oscillations found in experiments. State and parameter values of a point on a NS cycle in (32) are given in Table 4. The normal form coefficient of the NS cycle is . Since it is negative, we expect stable tori nearby.

Table 4: State and parameter values of a point on an NS- cycle in the Steinmetz-Larter model.
Variable | Value | Parameter | Value | Parameter | Value |
| 1.8609653 |  | 0.1631021 |  | 1.104 |
| 25.678306 |  | 1250 |  | 0.001 |
| 0.010838258 |  | 0.046875 |  | 0.71643356 |
| 0.094707061 |  | 20 |  | 0.1175 |
|  |  |  |  | 0.5 |

To start a time integration from this NS point, with a slightly supercritical parameter value, namely , we can use the commands

    X0=[1.8609653;25.678306;0.010838258;0.094707061];
    P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
    P0(7)=0.7167;
    hls=feval(@SteLar);
    options=odeset('RelTol',1e-8);
    [t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,1.104,
    0.001,0.7167,0.1175,0.5);
    size(t)
    size(y)
    plot(t(100000:339000),y(100000:339000,2));
    'press any key'
    pause
    plot(y(100000:339000,3),y(100000:339000,2));

After a transient the time-series exhibits modulated oscillations with two frequencies near the original limit cycle (see Fig. 7). This is a motion on a stable two-dimensional torus that arises from the Neimark-Sacker bifurcation. The commands can be executed by running the testrun `testtimeinteg.m`. The total number of computed points is 339161.

Figure 7: Movement on a stable torus in a 4-dimensional space as a function of time in the direction and as projected on the phase plane.

The testrun `testtimeintegJacobian.m`

    X0=[1.8609653;25.678306;0.010838258;0.094707061];
    P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
    P0(7)=0.7167;
    hls=feval(@SteLar);
    options=odeset('RelTol',1e-8);
    options = odeset('Jacobian',hls(3),'JacobianP',hls(4),
    'Hessians',hls(5),'HessiansP',hls(6));
    [t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,
    1.104,0.001,0.7167,0.1175,0.5);
    size(t)
    size(y)
    plot(t(100000:325000),y(100000:325000,2));
    'press any key'
    pause
    plot(y(100000:325000,3),y(100000:325000,2));

is a small modification of `testtimeinteg.m` in which the solver has access to the Jacobian matrix. The graphical output is similar, but the number of computed points is now 328197.