# square wave bursting
p v1=-1.2,v2=18,v3=12.,v4=17.4,gca=4.,gk=8.0,gl=2.0,vk=-84,vl=-60
p vca=120.0,I=45,gkca=.25,phi=.23,eps=.005,mu=.02,k=1
minf(v)=.5*(1.+tanh((v-v1)/v2))
winf(v)=.5*(1.+tanh((v-v3)/v4))
tauw(v)=1./cosh((v-v3)/(2*v4))
z(c)=c/(1.0+c)
v'= -(gca*minf(v)*(v-vca)+gk*w*(v-vk)+gl*(v-vl)+gkca*z(c)*(v-vk))+I
w'= 20*phi*(winf(v)-w)/tauw(v)
c'= 20*eps*(-mu*gca*minf(v)*(v-vca)-k*c)
i  v=-30,w=.0,c=1

aux axz=z(c)

@ dt=.05,total=200,meth=qualrk,xp=t,yp=v,xlo=0,xhi=200,ylo=-40,yhi=30

done