Examples of solving ODEs and visualizing their solutions in Mathematica 9.

This is a harmonic oscillator with cubic damping nonlinearity.  There are two functions of t, x[t], and y[t] and two first-order ODEs.  The initial condition is (x[0],y[0]) = (0.4,1.0).

soln=NDSolve[{x'[t]==-3y[t],y'[t]==x[t]-((y[t])^3),x[0]==0.4,y[0]==1.0},{x,y},{t,0,50}];

The output, which I called soln, is a rule that assigns to  x  and  y  the values of a function that interpolates the solution of the ODE at any values of  t  between 0 and 50.

soln

The functions  x  or  y   or their derivatives are obtained by using “Evaluate”.  Here are the values of  x’  at four values of  t:

Evaluate[x'[{2,3,5,40}]/.soln]

Here is a phase-plane plot of the solution.  It takes a looong time to converge to the origin like the neutral fixed points in 2030 class.

ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,50}]

You can plot  y  vs.  t, as long as t remains between 0 and 50, which is the domain of soln.

Plot[Evaluate[y[t]/.soln],{t,0,30}]

Here is the famous Van der Pol oscillator--a harmonic oscillator with nonlinear damping for large amplitude and anti-damping for small amplitude.  The damping is in the first-derivative term, with the parameter μ.  Try different values of  μ and different initial conditions.

μ=0.7;
VanDerPol=NDSolve[{y''[t]-μ(1-((y[t])^2))y'[t]+y[t]==0,y[0]==0.1,y'[0]==0.2},y,{t,0,60}];

The phase portrait is a parametric plot of  y  and  y’.

ParametricPlot[Evaluate[{y[t],y'[t]}/.VanDerPol],{t,0,50},Frame->True]

You can include the time axis in the plot.

ParametricPlot3D[Evaluate[{y[t],y'[t],t}/.VanDerPol],{t,0,40},AspectRatio->1]

The Duffing oscillator is a forced, damped, nonlinear oscillator that exhibits chaos in certain parameter regimes!

δ=0.2;β=1.0;α=1.0;γ=0.3;ω=0.5;
Duffing=NDSolve[{y''[t]+δ y'[t]+β y[t]+α ((y[t])^3)==γ Cos[ω t],y[0]==0,y'[0]==0},y,{t,0,40}];

Here’s a phase-plane plot.  To uncross the crossings, include the t axis using ParametricPlot3D.

ParametricPlot[Evaluate[{y[t],y'[t]}/.Duffing],{t,0,40},AspectRatio->1]

Here is the famous Lorenz chaotic “strange attractor” with typical parameters.

σ=10;ρ=28;β=8/3;
Lorenz=NDSolve[{x'[t]==σ(y[t]-x[t]),y'[t]==x[t](ρ-z[t])-y[t],z'[t]==x[t]y[t]-β z[t],x[0]==3,y[0]==15,z[0]==1},{x,y,z},{t,0,100},AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->40,MaxSteps->Infinity];

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.Lorenz],{t,0,100},PlotRange->All]

Here is a plot of  y  vs.  t.  You can see how often and how chaotically the trajectory passes from one wing of the butterfly to the other.

Plot[Evaluate[y[t]/.Lorenz],{t,0,100}]

Here’s a plot of   x  and  y  against  t.

ParametricPlot3D[Evaluate[{x[t],y[t],t}/.Lorenz],{t,0,100}]

Here is a mathematical model of spiking.  The differential equations are  x’(t) = f(x,y)  and  y’(t) = g(x,y),  where  f  and  g  are defined to allow the kind of bifurcations one sees in the example of cortical cells.  The parameter  z  controls the crossings between the two isoclines that produce two fixed points that undergo a tangent bifurcation.

(* f[x_,y_]:=a(y-1+c/(1+(x^2))); *)
f[x_,y_]:=a(y-0.5((x^3)-x+z));
(* g[x_,y_]:=-b (x +1)y/((1+(((x +1))^2))(0.2+9(y^2)));
g[x_,y_]:=-b (x +0.4)y/((1+(((x +0.4))^2))(0.2+9(y^2))); *)
g[x_,y_]:=-b (y-0.5(1-Tanh[10(x+0.4)]))/((0.3+(((x-0.5))^2))(1+2(((y-0.2))^2)));

Jac[x_,y_]:=Module[{s,t},{ D[f[s,t],{{s,t}}],D[g[s,t],{{s,t}}] }/.{s->x,t->y}];

FixedPoints:=Solve[{f[x,y]==0,g[x,y]==0},{x,y}];  (* doesn't finish *)

a=4;b=1;
(* z = 0.39 tangent bifurcation *)
z=0.35;
x0=-0.4;y0=0.3;
(* z=0.41;
x0=0.7;y0=0.01; *)
time=50;
soln=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,time}];
VectorField=VectorPlot[{f[x,y],g[x,y]},{x,-1.3,1.3},{y,-0.3,0.7}];
NullCurves=ContourPlot[{f[x,y]==0,g[x,y]==0},{x,-1.5,1.5},{y,-0.5,1},ContourStyle->{Darker[Red],Darker[Green]}];
Show[{VectorField,NullCurves},PlotRange->{{-1.3,1.3},{-0.2,0.6}}]

a=4;b=1;
(* z = 0.39 tangent bifurcation *)
z=0.28;
x0=-0.4;y0=0.3;
z=0.41;
x0=0.7;y0=0.01; 
time=50;
soln=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,time}];
VectorField=VectorPlot[{f[x,y],g[x,y]},{x,-1.3,1.3},{y,-0.3,0.7}];
NullCurves=ContourPlot[{f[x,y]==0,g[x,y]==0},{x,-1.5,1.5},{y,-0.5,1},ContourStyle->{Darker[Red],Darker[Green]}];
Trajectory=ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,time}];Show[{VectorField,NullCurves,Trajectory},PlotRange->{{-1.3,1.3},{-0.2,0.6}}]

Trajectory=Plot[Evaluate[x[t]/.soln],{t,0,time},PlotRange->All];Zvalue=Plot[z,{t,0,time},PlotRange->All,PlotStyle->Red];
Show[Trajectory,Zvalue,AspectRatio->0.2]

To simulate bursting, we couple the parameter  z  dynamically to  x  and  y.  It is a slowly varying quantity because of the small value of  ϵ  in its dynamical equation.  This contrived model causes  z  to increase if  x  is too big (at the upper fixed point) or decrease if  x  is too small (traversing the limit cycle).

When  z  hits its upper critical value where the two isoclines are tangent to each other,  x  suddenly drops into the limit cycle and  z  starts to decrease.  When  z  hits its lower critical value where the limit cycle bumps into the unstable fixed point,  x  suddenly shoots up to the stable fixed point and  z  starts to increase again.  Plotting the value of  z  shows how the system switches between spiking and quiescence at the two extrema of  z.

The parameter  w  controls the relative lengths of the spiking and quiescent states.  Increasing  w  lengthens the quiescent state.

f3[x_,y_,z_]:=a(y-0.5((x^3)-x+z));
(* g[x_,y_]:=-b (x +1)y/((1+(((x +1))^2))(0.2+9(y^2)));
g[x_,y_]:=-b (x +0.4)y/((1+(((x +0.4))^2))(0.2+9(y^2))); *)
g[x_,y_]:=-b (y-0.5(1-Tanh[10(x+0.4)]))/((0.3+(((x-0.5))^2))(1+2(((y-0.2))^2)));
h[x_,z_,w_]:=x-z-w;  (* h gets multiplied by ϵ in the differential equation. *)

Clear[z];
a=4;b=1;w=0.03;ϵ=0.007;  (* Play with a, b, w, and ϵ *)
x0=-0.4;y0=0.4;z0=0.4;
time=1000;
soln=NDSolve[{x'[t]==f3[x[t],y[t],z[t]],y'
[t]==g[x[t],y[t]],z'[t]==ϵ h[x[t],z[t],w],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,time},MaxSteps->Infinity];
(*Trajectory=ParametricPlot3D[Evaluate[{x[t],y[t],t}/.soln],{t,0,500}];*)
Trajectory=Plot[Evaluate[x[t]/.soln],{t,0,time},PlotRange->All];Zvalue=Plot[Evaluate[z[t]/.soln],{t,0,time},PlotRange->All,PlotStyle->Red];
Show[Trajectory,Zvalue,AspectRatio->0.2]

To simulate modulated bursting, we couple the parameter  w  dynamically to the others, in an ultra-slow way by scaling its dynamical equation by a very small number  δ  (much smaller than  ϵ).  This simple example causes  w  to decrease when the  x  value is greater than  xav  and to causes  w  to increase when the  x  value is less than  xav.  So the system approaches a bursting state in which the time-average of  x  is  xav,  which dictates a certain balance between the spiking and the quiescent states.

Plotting the value of  z  shows how the system switches between spiking and quiescence at the two extrema of  z.
Plotting the value of  w  shows how the average  x  value approaches  xav  as the bursting settles into a steady pattern.

Clear[w];
a=4;b=1;xav=0;ϵ=0.003;δ=0.00011;
x0=-0.4;y0=0.4;z0=0.4;w0=0.2;
time=5000;
soln=NDSolve[{x'[t]==f3[x[t],y[t],z[t]],y'[t]==g[x[t],y[t]],z'[t]==ϵ(x[t]-z[t]-w[t]),w'[t]==-δ(x[t]-xav),x[0]==x0,y[0]==y0,z[0]==z0,w[0]==w0},{x,y,z,w},{t,0,time},MaxSteps->Infinity];
(*Trajectory=ParametricPlot3D[Evaluate[{x[t],y[t],t}/.soln],{t,0,500}];*)
Trajectory=Plot[Evaluate[x[t]/.soln],{t,0,time},PlotRange->All];
Zvalue=Plot[Evaluate[z[t]/.soln],{t,0,time},PlotRange->All,PlotStyle->Red];
Wvalue=Plot[Evaluate[w[t]/.soln],{t,0,time},PlotRange->All,PlotStyle->Green];
Show[Trajectory,Zvalue,Wvalue,AspectRatio->0.2]

The potassium M-current neurobiological model

minf[V_]:=(1+Exp[(Vhalfm-V)/km])^(-1);
ninf[V_]:=(1+Exp[(Vhalfn-V)/kn])^(-1);
nMinf[V_]:=(1+Exp[(VhalfnM-V)/knM])^(-1);
tau[V_]:=Cbase+Camp Exp[-(((Vmax-V))^2)/(σ^2)];
tauM[V_]:=CMbase+CMamp Exp[-(((VMmax-V))^2)/((σM)^2)];
fV[V_,n_,nM_,Inp_]:= (Inp-gL(V-EL)-gNa minf[V](V-ENa)-gK n(V-EK)-gM*nM(V-EKM))/Current;
fn[V_,n_]:=(ninf[V]-n)/tau[V];
fnM[V_,nM_]:=ϵ(nMinf[V]-nM)/tauM[V];

Current=500;
Vhalfm=-40;Vhalfn=-40;km=9;kn=3.5;
Cbase=5;Camp=25;Vmax=-60;σ=30;
CMbase=100; CMamp=300; VMmax=-60; σM=25;
VhalfnM=-40; knM=10;
gL=0.3;gNa=15;gK=25;gM=42;
EL=10;ENa=110;EK=-100;EKM=-100;
(* The constants km and kn control the steepness of the positive-slope part of the zero sets of f and g. *)

By fixing  nM,  one obtains either a spiking state or a quiescent state.

nM=0.1;Inp=20;
V0=-40;n0=0.6;
tmax=1000;
soln=NDSolve[{V'[t]==fV[V[t],n[t],nM,Inp],n'[t]==fn[V[t],n[t]],V[0]==V0,n[0]==n0},{V,n},{t,0,tmax}];
VectorField=VectorPlot[{fV[V,n,nM,Inp],fn[V,n]},{V,-100,0},{n,-0.2,1.0}];
NullCurves=ContourPlot[{fV[V,n,nM,Inp]==0,fn[V,n]==0},{V,-100,0},{n,-0.2,1.0},ContourStyle->{Darker[Red],Darker[Green]}];
Trajectory=ParametricPlot[Evaluate[{V[t],n[t]}/.soln],{t,0,tmax}];Show[{VectorField,NullCurves,Trajectory},PlotRange->{{-100,0},{-0.2,1.0}}]

Now the slow variable  nM  is included dynamically to produce bursting.
I scaled all the equations in the ODE by 100 so that the time interval of computation could be reduced by a factor of 100.

Inp=20;CMbase=15000;ϵ=0.5;
Clear[nM];
V0=-40;n0=0.6;nM0=0.1;
tmax=1000;
soln=NDSolve[{V'[t]==100fV[V[t],n[t],nM[t],Inp],n'[t]==100fn[V[t],n[t]],nM'[t]==100fnM[V[t],nM[t]],V[0]==V0,n[0]==n0,nM[0]==nM0},{V,n,nM},{t,0,tmax},MaxSteps->Infinity];
Trajectory=Plot[Evaluate[V[t]/.soln],{t,0,tmax},PlotRange->All];nMvalue=Plot[100Evaluate[nM[t]/.soln],{t,0,tmax},PlotRange->All,PlotStyle->Red];
Show[Trajectory,AspectRatio->0.2]

qq[x_]:=-10x;
CC=NDSolve[{-y''[x]+qq[x]y[x]==0,y[0]==1,y'[0]==0},y,{x,0,1}];
SS=NDSolve[{-y''[x]+qq[x]y[x]==0,y[0]==0,y'[0]==1},y,{x,0,1}];
MM=Transpose[{Evaluate[{y[1],-y'[1]}/.CC][[1]],Evaluate[{y[1],-y'[1]}/.SS][[1]]}];
DD=DiagonalMatrix[Eigenvalues[MM]];
VV=Transpose[Eigenvectors[MM]];
cc[x_]:=Evaluate[y[x]/.CC];
ss[x_]:=Evaluate[y[x]/.SS];
u1[x_]:=VV[[1,1]]cc[x]+VV[[2,1]]ss[x];
u2[x_]:=VV[[1,2]]cc[x]+VV[[2,2]]ss[x];

Plot[{u1[x],u2[x],u1[x]u2[x]},{x,0,1}]

NIntegrate[u1[x]u2[x],{x,0,1}]
{Det[MM],MM.VV-VV.DD}

f1[x_,y_]:=-x(y+(x^2)-2x-1);
(* g1[x_,y_]:=y(x-0.8-0.1y); *)
g1[x_,y_]:=y(x-1);

x0=0.3;y0=0.6;time=100;
x1=-0.1;x2=2.2;y1=-0.1;y2=3;
NullCurves=ContourPlot[{f1[x,y]==0,g1[x,y]==0},{x,x1,x2},{y,y1,y2},ContourStyle->{Darker[Red],Darker[Green]}];VectorField=VectorPlot[{f1[x,y],g1[x,y]},{x,x1,x2+0.2},{y,y1,y2+0.25}];
soln=NDSolve[{x'[t]==f1[x[t],y[t]],y'[t]==g1[x[t],y[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,time}];
Trajectory=ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,time}];Show[{VectorField,NullCurves,Trajectory},PlotRange->{{x1,x2},{y1,y2}}]