Modified from Ch8.mws from Herod's web site by R.F.MurphyChapter 8, Sections 8.5 and 8.6An Introduction to the Mathematics of Biology, with Computer Algebra Modelsby Yeargers, Shonkwiler, & Herod. The Syntax is written for Maple 6Page 257Ena:=55: Ek:=-82: El:= -59: gkbar:=24.34: gnabar:=70.7:
gl:=0.3: vrest:=-69: cm:=0.001:alphan:=v-> 0.01*(10-(v-vrest))/(exp(0.1*(10-(v-vrest)))-1):
betan:=v-> 0.125*exp(-(v-vrest)/80):
alpham:=v-> 0.1*(25-(v-vrest))/(exp(0.1*(25-(v-vrest)))-1):
betam:=v-> 4*exp(-(v-vrest)/18):
alphah:=v->0.07*exp(-0.05*(v-vrest)):
betah:=v->1/(exp(0.1*(30-(v-vrest)))+1):pulse:=t->-20*(Heaviside(t-.001)-Heaviside(t-.002)):rhsV:=(t,V,n,m,h)->-(gnabar*m^3*h*(V-Ena) +
gkbar*n^4*(V-Ek) + gl*(V-El)+ pulse(t))/cm:
rhsn:=(t,V,n,m,h)-> 1000*(alphan(V)*(1-n) - betan(V)*n):
rhsm:=(t,V,n,m,h)-> 1000*(alpham(V)*(1-m) - betam(V)*m):
rhsh:=(t,V,n,m,h)-> 1000*(alphah(V)*(1-h) - betah(V)*h):inits:=V(0)=vrest,n(0)=0.315,m(0)=0.042, h(0)=0.608;sol:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric, output=listprocedure);Vs:=subs(sol,V(t));plot(Vs,0..0.02);sol20:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric);with(plots):J:=odeplot(sol20,[V(t),n(t)],0..0.02):display({J});pulse:=t->-2*(Heaviside(t-.001)-Heaviside(t-.002)):rhsV:=(t,V,n,m,h)->-(gnabar*m^3*h*(V-Ena) +
gkbar*n^4*(V-Ek) + gl*(V-El)+ pulse(t))/cm:sol2:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric);K:=odeplot(sol2,[V(t),n(t)],0..0.02,color=green):display({J,K});odeplot(sol20,[V(t),h(t)],0..0.02);odeplot(sol20,[m(t),h(t)],0..0.02);Figure 8.6.2a:=0.7; b:=0.8; c:=0.08;rhsx:=(t,x,y)->x-x^3/3-y;
rhsy:=(t,x,y)->c*(x+a-b*y);sol2:=dsolve({diff(x(t),t)=rhsx(t,x(t),y(t)),
diff(y(t),t)=rhsy(t,x(t),y(t)),x(0)=0,y(0)=-1},
{x(t),y(t)},type=numeric, output=listprocedure);xs:=subs(sol2,x(t)); ys:=subs(sol2,y(t));K:=plot([xs,ys,0..200],x=-3..3,y=-2..2,color=blue):
J:=plot({[V,(V+a)/b,V=-2.5..1.5],[V,V-V^3/3,V=-2.5..2.2]}):plots[display]({J,K});