# Traub spikes
itrb(v,m,h,n)=gna*h*m^3*(v-ena)+gk*n^4*(v-ek)+gl*(v-el)
# 
v'=-itrb(v,m,h,n) + i(t)
m'=am(v)*(1-m)-bm(v)*m
h'=ah(v)*(1-h)-bh(v)*h
n'=an(v)*(1-n)-bn(v)*n
# init data
init v=-66.6,m=.016,h=.995,n=.04

i(t)=i0+ip*heav(t-ton)*heav(toff-t)
am(v)=.32*(54+v)/(1-exp(-(v+54)/4))
bm(v)=.28*(v+27)/(exp((v+27)/5)-1)
ah(v)=.128*exp(-(50+v)/18)
bh(v)=4/(1+exp(-(v+27)/5))
an(v)=.032*(v+52)/(1-exp(-(v+52)/5))
bn(v)=.5*exp(-(57+v)/40)
#
# parameters
#
par i0=0,ip=2,ton=20,toff=50
par ek=-100,ena=50,el=-67
par gl=.1,gk=80,gna=100
@ total=200,dt=.1,meth=qualrk,toler=1e-5,atoler=1e-4
@ xhi=200,ylo=-100,yhi=50
done






x