restart; with(plots): with(DEtools): eq1 := C*diff(V(t),t)=I1-g[L]*(V(t)-E[L])-g[Na]*m[infinity](V(t))*(V(t)-E[Na])-g[K]*n(t)*(V(t)-E[K]); eq2 := diff(n(t),t)= (n[infinity](V(t))-n(t))/tau(V(t)); m1_inf := eval(subs({V[1/2]=-20,k=15},1/(1+exp((V[1/2]-V(t))/k)))); n1_inf := eval(subs({V[1/2]=-25,k=5},1/(1+exp((V[1/2]-V(t))/k)))); params := {C=1,I1=4.51,E[L]=-80,g[L]=8,g[Na]=20,g[K]=10,tau(V(t))=1,E[Na]=60,E[K]=-90, m[infinity](V(t))=m1_inf,n[infinity](V(t))=n1_inf}; eqn1a := eval(subs(params,eq1),{V(t)=V1,n(t)=n1}); eqn2a := rhs(eval(subs(params,eq2),{V(t)=V1,n(t)=n1}))=0; eqn1a := rhs(eval(subs(params,eq1),{V(t)=V1,n(t)=n1})); eqn2a := rhs(eval(subs(params,eq2),{V(t)=V1,n(t)=n1})); n1a := solve(eqn1a,n1); n2a := solve(eqn2a,n1); plot([n1a,n2a],V1=-89..0,n1=-0.1..0.5,axes=boxed); eqn1b := eval(subs(params,eq1)); eqn2b :=eval(subs(params,eq2)); dfieldplot([eqn1b,eqn2b],[V(t),n(t)],t=-5..5,V=-89..0,n=-0.1..0.5,arrows=small); plot1 := plot([n1a,n2a],V1=-89..0,n1=-0.1..0.5,axes=boxed,color=[blue,"DarkGreen"],thickness=2): plot2 := dfieldplot([eqn1b,eqn2b],[V(t),n(t)],t=-5..5,V=-89..0,n=-0.1..0.5,arrows=small): display(plot1,plot2); plot3 := DEplot([eqn1b,eqn2b],[V(t),n(t)],t=0..2,[[V(0)=-70,n(0)=0],[V(0)=-60,n(0)=0],[V(0)=-50,n(0)=0],[V(0)=-60,n(0)=0.1],[V(0)=-60,n(0)=-0.1]],V=-89..0,n=-0.1..0.5,arrows=small,linecolor=[black,blue,green,orange,cyan]); display(plot1,plot3);

restart; with(plots): with(DEtools): eq1 := C*diff(V(t),t)=I1-g[L]*(V(t)-E[L])-g[Na]*m[infinity](V(t))*(V(t)-E[Na])-g[K]*n(t)*(V(t)-E[K]); eq2 := diff(n(t),t)= (n[infinity](V(t))-n(t))/tau(V(t)); m1_inf := eval(subs({V[1/2]=-20,k=15},1/(1+exp((V[1/2]-V(t))/k)))); n1_inf := eval(subs({V[1/2]=-25,k=5},1/(1+exp((V[1/2]-V(t))/k)))); params1 := {C=1,E[L]=-80,g[L]=8,g[Na]=20,g[K]=10,tau(V(t))=1,E[Na]=60,E[K]=-90, m[infinity](V(t))=m1_inf,n[infinity](V(t))=n1_inf}; eqn1c := rhs(eval(subs(params1,eq1),{V(t)=V1,n(t)=n1}))=0; eqn2c := rhs(eval(subs(params1,eq2),{V(t)=V1,n(t)=n1}))=0; n1c := solve(eqn1c,n1); n2c := solve(eqn2c,n1); eqn1d := eval(subs(params1,eq1)); eqn2d :=eval(subs(params1,eq2)); animate(plot,[[n1c,n2c],V1=-89..0,n1=-0.1..0.5,axes=boxed],I1=-20..20); eqn1b := eval(subs(params1,eq1)); eqn2b :=eval(subs(params1,eq2)); animate(DEplot,[(eval([eqn1d,eqn2d],I1=I2),[V(t),n(t)],t=0..2,[[V(0)=-70,n(0)=0],[V(0)=-60,n(0)=0],[V(0)=-50,n(0)=0],[V(0)=-60,n(0)=0.1],[V(0)=-60,n(0)=-0.1]],V=-89..0,n=-0.1..0.5,arrows=small,linecolor=[black,blue,green,orange,cyan])],I2=-40..20); plotseq1 := []: for I1 from -30 by 1 to 20 do plot1d := plot([n1c,n2c],V1=-89..0,n1=-0.1..0.5,axes=boxed,color=[blue,"DarkGreen"],thickness=1): plot2d := DEplot([eqn1d,eqn2d],[V(t),n(t)],t=0..5, [[V(0)=-80,n(0)=0],[V(0)=-70,n(0)=0],[V(0)=-60,n(0)=0],[V(0)=-50,n(0)=0],[V(0)=-40,n(0)=0], [V(0)=-80,n(0)=0.1],[V(0)=-70,n(0)=0.1],[V(0)=-60,n(0)=0.1],[V(0)=-50,n(0)=0.1],[V(0)=-40,n(0)=0.1], [V(0)=-80,n(0)=-.1],[V(0)=-70,n(0)=-.1],[V(0)=-60,n(0)=-.1],[V(0)=-50,n(0)=-0.1],[V(0)=-40,n(0)=-0.1]], V=-89..0,n=-0.1..0.5,arrows=small,linecolor=[black,orange,green,cyan,violet,black,orange,green,cyan,violet,black,orange,green,cyan,violet]): plotseq1 := [op(plotseq1),display(plot1d,plot2d)]; end do: display(plotseq1,insequence=true);