# bistability in autocatalytic rxn
param dt=0.0001 kf1=0.001 kb1=10 kf2=250
d' = -(f1-b1)/dt
dp2' = (f1-b1)/dt
p' = -2*(f1-b1)/dt + (f2-f3)/dt
init d=1 dp2=0 p=100
markov f1 2
{0} {kf1*d*p*(p-1)}
{1/dt-kf1*d*p*(p-1)} {0}
markov b1 2
{0} {kb1*dp2}
{1/dt-kb1*dp2} {0}
markov f2 2
{0} {kf2*dp2}
{1/dt-kf2*dp2} {0}
markov f3 2
{0} {p}
{1/dt-p} {0}
@ meth=euler bound=1000 maxstor=10000000 dt=0.0001
done