plotting – 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:

enter code here

How do I modify the above code, so that it solves the below system?
enter image description here

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!