# bruss_gil.ode
# gillespie algorithm for brusselator
#
# here are the four reactions
# A  -> U (c1)
# B+U -> V (c2)
# 2 U + V -> 3 U (c3)
# U -> * (c4)
# 
par c1a=5000,c2b=50,c3=.00005,c4=5
init u=1000,v=2000
#  compute the cumulative reaction rates
# 
p1=c1a
p2=p1+c2b*u
p3=p2+c3*u*v*(u-1)/2
p4=p3+c4*u
# choose random # to determine which reaction
# note that p4 is the total of all the rates 
s2=ran(1)*p4
# z1,z2,z3,z4 are either 0 or 1 according to whether the reaction
# took place or not
# if s2<p1, then z1=1; if p1<=s2<p2 then reaction 2 takes place,...
#
z1=(s2<p1)
z2=(s2<p2)&(s2>=p1)
z3=(s2<p3)&(s2>=p2)
# finally, if s2 >= p3, then reaction 4
z4=(s2>=p3)
# time for next reaction
tr'=tr-log(ran(1))/p4
u'=max(0,u+z1-z2+z3-z4)
v'=max(0,v+z2-z3)
# this plots a phaseplane of the variables - it does 1000000 reactions
# and ploats out every thousandth
@ bound=100000000,meth=discrete,total=1000000,nout=1000
@ xp=u,yp=v
@ xlo=0,ylo=0,xhi=10000,yhi=10000
done