#mackey2.ode
# this is the kernel:
xi(t)=heav(t-1)*heav(tmax-t)/(tmax-1)
# I have changed theta a bit to make sure it isnt undefined at t=0
# 
th(t)=1/(t+eps)^3
par eps=.001
#
g(y)=y/(1+y^n)
#  max(u,0) == heav(u)*u
#
v'=gamma*(e-v)-beta0*g(f0*int{max(v-th(t-t'),0)*xi(t-t')})
p gamma=10,n=3,f0=8,beta0=600,tmax=1.625,e=1.6
init v=.9
@ TOTAL=20,DT=0.01,XLO=0,XHI=20,YLO=-30,YHI=10,VMAXPTS=400
done