#  hhh.ode
init v=-65,  m=.05, h=0.6  n=.3
par iapp=0
par vna=50  vk=-77  vl=-54.4  gna=120  gk=36  gl=0.3 
am(v)=.1*(v+40)/(1-exp(-(v+40)/10))
bm(v)=4*exp(-(v+65)/18)
ah(v)=.07*exp(-(v+65)/20)
bh(v)=1/(1+exp(-(v+35)/10))
an(v)=.01*(v+55)/(1-exp(-(v+55)/10))
bn(v)=.125*exp(-(v+65)/80)
minf(v)=am(v)/(am(v)+bm(v))
taum(v)=1/(am(v)+bm(v))
hinf(v)=ah(v)/(ah(v)+bh(v))
tauh(v)=1/(ah(v)+bh(v))
ninf(v)=an(v)/(an(v)+bn(v))
taun(v)=1/(an(v)+bn(v))
v'= - gna*h*(v-vna)*m^3-gk*(v-vk)*n^4-gl*(v-vl)+iapp
m'=(minf(v)-m)/taum(v)
h'=(hinf(v)-h)/tauh(v)
n'=(ninf(v)-n)/taun(v)

# track the currents
aux ina=gna*(v-vna)*h*m^3
aux ik=gk*(v-vk)*n^4
aux il=gl*(v-vl)

@ dt=.05,total=100,meth=qualrk,xp=t,yp=v,xlo=0,xhi=100,ylo=-80,yhi=20.,bound=500001,maxstor=5000001

done