# Solve a system of three stiff ODEs and plot graphs of it - Oregonator

I want to plot a system of three stiff ODEs (Oregonator model). That model describes a chemical oscillator. I do not have much experience with plotting ODEs, but I have obtained a Mathematica file that does almost the same as I want, but only with a simplified version of the Oregonator. That simplified version has two ODEs, while my system has three ODEs.

I have the following code:

ε=4*10^-2;
δ=4*10^-4;
q=8*10^-4;
f=1;

{xsol,zsol}=NDSolve[{εx'[t]==x[t](1-x[t])+(f(q-x[t])z[t])/(q+x[t]),z'[t]==x[t]-z[t],
x[0]==.00012,z[0]==.00576},{x,z},{t,0,40},MaxSteps->Infinity]


This system describes the following model:

$$\epsilon \frac{dx}{dt}=x(1-x)+\frac{f(q-z)}{q+x}z$$ $$\delta \frac{dz}{dt}=x-z$$

How do I modify the above code, so that it solves the below system?

I plot the solutions like this:

Plot[Evaluate[x[t]/.xsol],{t,0,40},PlotRange->All,PlotStyle->{Thick,Blue}]


I already tried adding ysol and y[t] and y[0], but that didn't work for me. Probably I did something wrong, but I'm really stuck. The initial value is y[0]=0.375.

How should I modify the code with the NDSolve? Thanks a lot!

REQUEST: it is appreciated if sb replaces the images of the system with the the correct MathJax-notation. I don't have experience with MathJax.

• Should the last term in the first equation of the 3-species system be $x(1-x)$ instead of $x(1-z)$? – Michael E2 May 23 '20 at 15:45
• @MichaelE2 That's indeed true. Sharp eye do you have ;). Thanks for being so attentive leaving a comment! I'll change the question and the answer so that it represents the correct system. – ralphjsmit May 25 '20 at 7:03

How do I modify the above code, so that it solves the below system?

You did not have IC for $$y(t)$$, so I made one up.

Clear["Global*"];
ε = 4*10^-2;
δ = 4*10^-4;
q = 8*10^-4;
f = 1;

ode1 = ε x'[t] == q y[t] - x[t]*y[t] + x[t] (1 - x[t]);
ode2 = δ y'[t] == -q y[t] - x[t]* y[t] + f*z[t];
ode3 = z'[t] == x[t] - z[t];

ic = {x[0] == .00012, y[0] == .375, z[0] == .00576};

{xSol, ySol, zSol} = NDSolveValue[{ode1, ode2, ode3, ic}, {x, y, z}, {t, 0, 40}]


To plot, one way

p[f_] := Plot[f, {t, 0, 40}, PlotRange -> All];
Grid[{{p[xSol[t]], p[ySol[t]], p[zSol[t]]}}, Frame -> All]


• Hi @Nasser thanks for also updating the images and of course the answer in the first place. In the first image for x[t], some parts of the graph are not displayed. Do you know how that comes and how I can display the first graph completely? So that I have an idea what the maximum is for x[t]. – ralphjsmit May 25 '20 at 9:54
• @ralphsmit updated. You can use PlotRange->All` without this option, Mathematica makes its own mind what is a "good" range. – Nasser May 25 '20 at 9:56
• @Nassar Awesome! Thanks for telling met that! – ralphjsmit May 25 '20 at 9:57