@ autoeval=0,total=400,nout=10
par vmin=0,vmax=5
b_(v)=a*v/(c+v)
par a=8,c=1,k=.75,u=1
table v_ % 101 0 100 ran(1)
table v2 % 101 0 100 sum(0,t)of(v_(i'))
table v % 101 0 100 vmax*v2(t)/v2(100)
table b % 101 0 100 b_(v(t))
ps0=0
ps[1..100]=ps[j-1]+y[j-1]
pbs0=b(0)*y0
pbs[1..100]=pbs[j-1]+b([j])*y[j]
x'=k-u*x-x*pbs100
y[0..100]'=y[j]*(b([j])*x-u-v([j])+s*b([j])*ps[j]-s*(pbs100-pbs[j]))
par s=2
init y[0..100]=.001
init x=1
aux pbs_100=pbs100
aux ps_100=ps100
done