# phasepul.ode
# pulsatile phase model
r(x)=-a*sin(x)
p(x)=exp(-beta*(1-cos(x)))
par a=3,beta=5
x1'=w1+p(x2)*r(x1)
x2'=w2+p(x1)*r(x2)
par w1=1,w2=3
@ total=100
@ fold=x1,fold=x2
@ xp=x1,yp=x2,xlo=0,ylo=0,xhi=6.3,yhi=6.3
done