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

Question

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.

Answer 1

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]