# here is gillespie of LV equations
# xb+y1 -> 2 Y1 (c1)
# y1 + y2 -> 2 Y2 (c2)
# y2 -> Z (c3)
#  
par c1x=10,c2=.01,c3=10
init tr=0,y1=1000,y2=1000
# compute cumulative sum of rates
a1=c1x*y1
a2=a1+y1*y2*c2
a3=a2+y2*c3
#  now choose the random number
q=ran(1)*a3
#  choose the reaction - zj=1 iff that reaction is run!
# choose 1 if q<a1, choose 2 if q < a2 and >= a1
# choose 3 if q>a2
z1=(q<a1)
z2=(q<a2)*(1-z1)
z3=(q>a2)
# time for next reaction
tr'=tr-log(ran(1))/a3
# make sure no one goes to 0
y1'=max(1,y1+z1-z2)
y2'=max(1,y2+z2-z3)
#      
@ bound=100000000,meth=discrete,total=1000000,njmp=1000
@ xp=y1,yp=y2
@ xlo=0,ylo=0,xhi=4000,yhi=4000
done