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]=-30,k=7},1/(1+exp((V[1/2]-V(t))/k)))); n1_inf := eval(subs({V[1/2]=-45,k=5},1/(1+exp((V[1/2]-V(t))/k)))); params1 := {C=1,E[L]=-78,g[L]=1,g[Na]=4,g[K]=4,tau(V(t))=1,E[Na]=60,E[K]=-90, m[infinity](V(t))=m1_inf,n[infinity](V(t))=n1_inf}; #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..20,n1=-0.1..1,axes=boxed],I1=43..49); eqn1b := eval(subs(params1,eq1)); eqn2b :=eval(subs(params1,eq2)); animate(DEplot,[(eval([eqn1d,eqn2d],I1=I2),[V(t),n(t)],t=0..8,[[V(0)=-70,n(0)=.1],[V(0)=-50,n(0)=0.3],[V(0)=-50,n(0)=0.5],[V(0)=-50,n(0)=0.7],[V(0)=-60,n(0)=-0.1]],V=-89..20,n=-0.1..1,arrows=small,linecolor=[black,blue,green,orange,cyan])],I2=43..49); plotseq1 := []: for I1 from 43 by 0.5 to 49 do plot1d := plot([n1c,n2c],V1=-89..20,n1=-0.1..1,axes=boxed,color=[blue,"DarkGreen"],thickness=1): plot2d := DEplot([eqn1d,eqn2d],[V(t),n(t)],t=0..5, [[V(0)=-80,n(0)=0.3],[V(0)=-70,n(0)=0.3],[V(0)=-60,n(0)=0.3],[V(0)=-50,n(0)=0.3],[V(0)=-40,n(0)=0.3], [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)=0.5],[V(0)=-70,n(0)=0.5],[V(0)=-60,n(0)=0.5],[V(0)=-50,n(0)=0.5],[V(0)=-40,n(0)=0.5]], V=-89..20,n=-0.1..1,stepsize=0.05,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);