Differential Equations

Math 2326

Helmut Knaust, Department of Mathematical Sciences, UTEP, El Paso TX 79968
hknaust@utep.edu

Last Edits:2/11/2026

Evaluate a command by pressing  ⇧ and  ⌤ at the same time. Evaluate the "Preparations" bracket first :

Preparations

SlopeField[f_,{iv_,st_,end_},{dv_,dst_,dend_}]:=Module[{Fct,t,y},Fct=Function[{iv,dv},f];VectorPlot[{1,Fct[t,y]},{t,st,end},{y,dst,dend},Axes->True,VectorStyle->Arrowheads[0],AxesLabel->{iv,dv},BaseStyle->{14,Bold},VectorScale->{Tiny,Small,None},PlotLabel->Row[{dv,"' = ",Fct[iv,dv]}]]]
Euler[rhs_,t0_,y0_,tend_,Δt_]:=Module[{G,f,eq,sol,pts,rhs2,min,max,opt},rhs2=rhs/.y->y[t];eq={y'[t]==rhs2,y[t0]==y0};f[a_,b_]:=rhs/.{t->a,y->b};
G[{t_,y_}]:={t+Δt,y+Δt f[t,y]/.{t->t,y->y}};sol=DSolve[eq,y,t];pts=NestList[G,{t0,y0},Ceiling[(tend-t0)/Δt]];min=Min[pts[[All,2]]];max=Max[pts[[All,2]]];z[t_]=Evaluate[y[t]/.sol[[1]]];
opt=AbsoluteOptions[Plot[z[t],{t,t0,tend},PlotRange->Automatic],PlotRange][[1,2,2]];
Plot[z[t],{t,t0,tend},PlotStyle->{Blue,AbsoluteThickness[2]},ImageSize->800,GridLines->Automatic,AxesLabel->{"t","y(t)"},PlotRange->{Min[opt[[1]],min]-1,Max[opt[[2]],max]+1},Epilog->{{PointSize[Medium],Red,Point[pts]},{Black,PointSize[Large],Point[{t0,y0}]}}]];
PhaseLine[f_,{t_,ts_,te_},{y_,ys_,ye_}]:=Module[{F},F=Function[y,f];VectorPlot[{1,F[y]}/Norm[{1,F[y]}], {t,ts,te}, {y, ys,ye},Axes->{True,True},LabelStyle->{Black,Bold,{FontSize->18}},PlotRange->{{ts,te},{ys,ye}},AxesOrigin->{0,0},AspectRatio->1/GoldenRatio,AxesLabel->{t,y},ImageSize->800,VectorScale->{Tiny,Small,None},Epilog->{Table[{If[F[y]>0,Red,Blue],PointSize[Medium],Point[{0,y}]},{y,ys,ye,(ye-ys)/500}],Select[{Green,AbsoluteThickness[2],Line[{{ts,y},{te,y}}]}/.Solve[F[y]==0],(ys≤#[[3,1,1,2]]≤ye)&]}]]
BifurcationDiagram[f_,{y_,ys_,ye_},{p_,ps_,pe_}]:=Module[{F},F=Function[{y,p},f];ContourPlot[F[dep,par],{par,ps,pe},{dep,ys,ye},Contours->{0},ContourStyle->Directive[White,AbsoluteThickness[4]],Axes->True,AxesLabel->{p,y},PlotRange->All,ContourShading->{Blue,Red},AspectRatio->1/GoldenRatio,ImageSize->800, PerformanceGoal->"Quality",LabelStyle->{Black,Bold,{FontSize->18}}]]
EulerSystem[rhs1_,rhs2_,t0_,x0_,y0_,tend_,Δt_]:=Module[{G,f,g,eq,sol,pts,pts2,Rhs1,Rhs2,p1,p2},f[a_,b_,c_]:=rhs1/.{t->a,x->b,y->c};g[a_,b_,c_]:=rhs2/.{t->a,x->b,y->c};
G[{t_,x_,y_}]:={t+Δt,x+Δt f[t,x,y]/.{t->t,x->x,y->y},y+Δt g[t,x,y]/.{t->t,x->x,y->y}};pts=NestList[G,{t0,x0,y0},Ceiling[(tend-t0)/Δt]];pts2=Transpose[Transpose[pts][[2;;3]]];p1=ListPlot[pts2,Joined->True,ImageSize->800,AxesLabel->{x,y},LabelStyle->{Black,Bold,{FontSize->18}},PlotStyle->{Thick,Blue},PlotRange->All,Epilog->{PointSize[Large],Point[{x0,y0}]}];
Rhs1=rhs1/.{x->x[t],y->y[t]};Rhs2=rhs2/.{x->x[t],y->y[t]};eq={x'[t]==Rhs1,y'[t]==Rhs2,x[t0]==x0,y[t0]==y0};sol=NDSolve[eq,{x,y},{t,t0,tend}];p2=ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,t0,tend},PlotStyle->{Lighter[Blue],Thin,Dashed}];Show[{p1,p2},PlotLabel->"Euler approximation and Mathematica approximation (dashed)"]];
Linearization[f_,g_,x_,y_,pt_]:=Module[{F,G,A,z},F=Function[{x,y},f];G=Function[{x,y},g];z=If[F[pt[[1]],pt[[2]]]==0&&G[pt[[1]],pt[[2]]]==0," "," NOT "];Print ["The point ",pt," is",z,"an equilibrium point.\nThe linearized system at ",pt," has the form Y'(t)= ",A={{D[F[x,y],x],D[F[x,y],y]},{D[G[x,y],x],D[G[x,y],y]}}/.{x->pt[[1]],y->pt[[2]]}," · Y(t).\nThe eigenvalues are \!\(\*SubscriptBox[\(λ\), \(1\(\\ \)\)]\) = ",Eigenvalues[N[A]][[1]]," and (Subscript[λ, 2]) = ", Eigenvalues[N[A]][[2]],"."]]
NullClinePlot[v1_,v2_,xr_,yr_]:=Module[{p1,p2,p3,s,F1,F2},F1=Function[{x,y},v1];F2=Function[{x,y},v2];
p1=ContourPlot[v1,xr,yr,PlotPoints->75,Contours->{0},PlotRange->All,ContourStyle->Directive[Red,Opacity[1],AbsoluteThickness[6]],ContourShading->{{Opacity[.4],GrayLevel[1]},{Opacity[.4],GrayLevel[0.1]}}];
p2=ContourPlot[v2,xr,yr,PlotPoints->75,Contours->{0},PlotRange->All,ContourStyle->Directive[Blue,Opacity[1],AbsoluteThickness[4]],ContourShading->{{Opacity[.4],GrayLevel[1]},{Opacity[.4],GrayLevel[0.2]}}];
p3=VectorPlot[If[{v1,v2}=={0,0},{0,0},{v1,v2}/Norm[{v1,v2}]], xr, yr,Axes->True,VectorScale->{Tiny,Small,None},BaseStyle->GrayLevel[.5]];
Show[p2,p1,p3,ImageSize->800,PlotLabel->StringForm["````\t\t````",Style[xr[[1]],Red,18,Bold],Style["'=0",Red,18,Bold],Style[yr[[1]],18,Bold,Blue],Style["'=0",18,Bold,Blue]],Axes->True,AxesLabel->{xr[[1]],yr[[1]]},LabelStyle->{Black,Bold,{FontSize->18}},Frame->False,PlotRangePadding->.01Min[xr[[3]]-xr[[2]],yr[[3]]-yr[[2]]],
Epilog->{Yellow,AbsolutePointSize[10],Point[s=Select[{x,y}/.Solve[{F1[x,y]==0,F2[x,y]==0},{x,y}],Element[#[[1]],Reals]&&Element[#[[2]],Reals]&]]}]]

Chapter 1

1.3 Slope Fields

The following command plots a slope field; you have to provide ranges for both the independent and dependent variable.

SlopeField[y^2 +Cos[t^2],{t,-2,2},{y,-3,3}]

1.4 Euler' s Method

The following command plots the exact solution to a differential equation (in blue) and its Euler approximation (in red).  The starting point is indicated in black. This only works if Mathematica can find an exact solution.

A typical example is 
		Euler[t Sin[t]/(y^2 + 10), 0, 5, 20, 0.5]
Inputs needed are 
	(i) the right-hand side of the differential equation (t Sin[t]/(y^2 + 10) in the example above)
	(ii and iii) the initial conditions (Subscript[t, 0]) and (Subscript[y, 0]) (in the example  0 and 5 respectively)
	(iv) the upper bound for t in the viewing window (in our case 20)
	(v) the step size for the Euler approximation (in the example 0.25).

Euler[t Sin[t]/(y^2+10),0,5,20,.5]

1.6 Phaseline Analysis for Autonomous Differential Equations

The following command plots the phaseline, the slope field and the equilibrium solutions of an autonomous differential equation.

A typical example is 
	Phaseline[.2 (z + .2) (z - 3.8)^2 (z - 2.1), {t, -2, 2},{z, -2, 4}]

Inputs needed are the right-hand side of the differential equation, the name of the independent variable and its start and endpoint, and the dependent variable name with start and endpoint.

Red (blue) portions of the phaseline indicate regions where the solutions are increasing (decreasing); equilibrium solutions are plotted in green.

PhaseLine[.2(z+.2)(z-3.8)(z-2.1),{t,-1,3},{z,-1,4}]

1.7 Bifurcation Diagrams

The command below produces a bifurcation diagram for a one-parameter family of autonomous differenial equations. 

Inputs needed are the right-hand side of the family of differential equations, the name of the dependent variable and its start and endpoint, and the parameter name with start and endpoint.
  
Red (blue) portions of the diagram indicate regions where the solutions are increasing (decreasing); the locus of the equilibrium points is plotted in white.

Clear[y,β];
BifurcationDiagram[(y+2β)(y^3-β),{y,-1,1},{β,-0.5,0.5}]

The animation below displays  a "Subcritical Pitchfork Bifurcation".

Manipulate[
ContourPlot[r*y-4*y^3+h,{r,-1,1},{y,-1,1},Contours->{0},ContourStyle->{White,AbsoluteThickness[4]},ContourShading->{Blue,Red},Axes->True,AxesLabel->{"r","y"},BaseStyle->{14,Bold},ImageSize->450,PlotLabel->Row[{"r y+4(y^3)",If[Chop[h]==0,"",PaddedForm[Chop[h],{3,2},NumberSigns->{"-","+"}]]}],PerformanceGoal->"Quality"],
{{h,0,""},-0.2,0.2,0.01}]

Chapter 2

2.1 A Solution Curve and its Projections

Clear[R,F];deq:={R'[t]==2R[t](1-R[t]/10)-R[t]  F[t],F'[t]==-F[t]+2R[t]  F[t],R[0]==1,F[0]==.2};Column[deq,ItemStyle->{{FontSize->16}}]
sol=NDSolve[deq,{R,F},{t,0,25}];
R=R/.sol[[1]]; F=F/.sol[[1]];

Manipulate[ParametricPlot3D[{R[t],k F[t], t},{t,0,15},AxesLabel->{"R","F",t},PlotStyle->{Thick,Blue},PlotRange->{{0,2.5},{0,6},{0,15}},BoxRatios->{1,1,1},FaceGrids->{{0,-1,0}},FaceGridsStyle->Directive[Gray,Thin],BaseStyle->{14,Bold}],{k,1,0,-.1}]

Manipulate[ParametricPlot3D[{k R[t], F[t], t},{t,0,15},AxesLabel->{"R","F",t},PlotStyle->{Thick,Blue},PlotRange->{{0,2.5},{0,6},{0,15}},BoxRatios->{1,1,1},FaceGrids->{{-1,0,0}},FaceGridsStyle->Directive[Gray,Thin],BaseStyle->{14,Bold}],{k,1,0,-.1}]

Manipulate[ParametricPlot3D[{R[t],F[t],k t},{t,0,15},AxesLabel->{"R","F",t},PlotStyle->{Thick,Blue},PlotRange->{{0,2.5},{0,6},{0,15}},BoxRatios->{1,1,1},BaseStyle->{14,Bold},FaceGrids->{{0,0,-1}},FaceGridsStyle->Directive[Red,Gray]],{k,1,0,-.1}]

ParametricPlot3D[{{R[t],0,t},{0, F[t], t},{R[t],F[t],0}},{t,0,15},AxesLabel->{"R","F",t},PlotStyle->{{Thick,Blue},{Thick,Green},{Thick,Red}},PlotRange->{{0,2.5},{0,6},{0,15}},BoxRatios->{1,1,1},BaseStyle->{14,Bold}]

2.2.1 The Glider Example

Here is the autonomous system of differential equations the British engineer F.W. Lanchester used in 1907 to model the motion of a glider. d denotes the drag coefficient of the glider.

glider:={v'[t]==-Sin[θ[t]]-d (v[t])^2,θ'[t]==(-Cos[θ[t]]+(v[t])^2)/v[t]};Column[glider,ItemStyle->{{FontSize->16}}]

We will first study the model without drag (in particular this means that the glider will stay in the air forever).

d=0;Column[glider,ItemStyle->{{FontSize->16}}]
s=Table[NDSolve[Union[glider,{v[0]==k,θ[0]==0}],{v,θ},{t,-10,10}],{k,1,2.8,.1}];
ParametricPlot[Evaluate[{v[t],θ[t]}/.s],{t,-5,5},AxesLabel->{v,θ},AspectRatio->1,PlotRange->{{0,3},{-π,π}},ImageSize->800, PlotStyle->AbsoluteThickness[2],LabelStyle->{Black,Bold,{FontSize->18}},Ticks->{Automatic,{-π,-π/2,0,π/2,π}}]

Here is the same model with a small drag coefficient :

d=0.1; Column[glider,ItemStyle->{{FontSize->16}}]
s=Table[NDSolve[Union[glider,{v[0]==k,θ[0]==0}],{v,θ},{t,0,20}],{k,1,5,.25}];
ParametricPlot[Evaluate[{v[t],θ[t]}/.s],{t,0,20},AxesLabel->{v,θ},AspectRatio->1,ImageSize->800, Ticks->{Automatic,Table[k π,{k,-1/2,4,1/2}]},PlotStyle->AbsoluteThickness[2],LabelStyle->{Black,Bold,{FontSize->18}},GridLines->{Automatic,Table[k π,{k,0,4,1}]},PlotRange->{{0,5},{-π/2,9π/2}}]

Here is  the Mathematica code to visualize flight paths :

d=0.1;Column[glider,ItemStyle->{{FontSize->16}}]
s=NDSolve[Union[glider,{v[0]==4,θ[0]==-π/4}],{v,θ},{t,0,30}];
Δt=.15;vθ=Flatten[Table[Evaluate[{v[t],θ[t]}/.s],{t,0,30,Δt}],1];
xy=Prepend[Drop[Map[({#[[1]] Δt Cos[#[[2]]],#[[1]]Δt Sin[#[[2]]] }&),vθ],1],{0,0}];
ListPlot[acc=Accumulate[xy],Joined->False,Axes->False,Frame->True,FrameTicks->None,FrameLabel->{"horizontal distance (x)","height above ground (y)"},AspectRatio->1/4,LabelStyle->{Black,Bold,{FontSize->14}},ImageSize->800,PlotStyle->{Blue,Thick}]

2.2.2 Bungee Jumping

The following second order non-linear differential equation models a bungee jump:

	m y''[t] + m g + b[y[t]] + β y'[t] == 0
	
y[t] denotes the position of the jumper at time t;  m is the mass of the jumper,  g is the gravitational constant; β is the damping coefficient. 
A bungee cord acts like a spring when stretched, but it has no restoring force when "compressed". The restoring force of the bungee cord b(y)  is thus modeled as follows: 
		b(y)=k y 	if y<0 	
	and 	b(y)=0		 if y≥0.

g=32;m=160/32;β=1;k=15;b[y_]=k y UnitStep[-y];
Plot[b[y],{y,-20,20},ImageSize->450,PlotStyle->AbsoluteThickness[3],Axes->{False,True},Frame->True,FrameLabel->{y,"b(y)"},LabelStyle->{Black,Bold,{FontSize->18}}]

Below are the graphs of the position of the jumper, the jumper's  velocity and the jumper's acceleration. The last graph shows the solution in the phase plane.

g=32;m=160/32;β=1;k=15;b[y_]=k y UnitStep[-y];
sol=NDSolve[{m y''[t]+ m g +b[y[t]]+β y'[t]==0,y[0]==100,y'[0]==0},y,{t,0,50},MaxStepFraction->1/100];
Plot[Evaluate[{y[t],y'[t]}/.sol],{t,0,10},PlotPoints->500,PlotRange->All,PlotStyle->{Directive[Red,AbsoluteThickness[3]],Directive[Blue,Opacity[1],AbsoluteThickness[1]]},ImageSize->600,PlotLabel->StringForm["``\t``",Style["y(t)",Red,18,Bold],Style["y'(t)",18,Bold,Blue]],LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{t,None}]
Plot[Evaluate[{y[t],y''[t]}/.sol],{t,0,10},PlotPoints->500,PlotRange->All,PlotStyle->{Directive[Red,AbsoluteThickness[3]],Directive[Orange,Opacity[1],AbsoluteThickness[1]]},ImageSize->600,PlotLabel->StringForm["``\t``",Style["y(t)",Red,18,Bold],Style["y''(t)",18,Bold,Orange]],LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{t,None}]
sol2=NDSolve[{y'[t]==v[t],v'[t]+ g +b[y[t]]/m+β v[t]/m==0,y[0]==100,v[0]==0},{y,v},{t,0,50},MaxStepFraction->1/100];ParametricPlot[Evaluate[{y[t],v[t]}/.sol2],{t,0,30},PlotPoints->500,PlotRange->All,PlotStyle->Directive[Black,AbsoluteThickness[3]],ImageSize->600,LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{y,v}]

2.4.1 Euler' s Method for Systems

The following command plots the Euler approximation of a solution in the phase plane.  The starting point is indicated in black. 

A typical example is 
		EulerSystem[Sin[x]/(y^2 + x),y^2-x, 0, 5,-2,20,0.25]
Inputs needed are 
	(i and ii) the right-hand sides of the differential equations 
	(iii-v) the initial conditions (Subscript[t, 0]) ,(Subscript[x, 0])and (Subscript[y, 0])  
	(vi) the upper bound for t in the viewing window (in our case 20)
	(vii) the step size for the Euler approximation (in the example 0.25).

EulerSystem[Sin[x]/(y^2+1),0.1y^2-x,0,5,1,20,0.35]

Below is a non - autonomous example. (Note that solutions of non-autonomous can self-intersect in the phase plane!)

EulerSystem[y -1,-x-.2y (t-5)+1,0,-1,0,11,0.2]

2.4.2 The Existence and Uniqueness Theorem Applied to Phase Portraits of Autonomous Systems

z[a_]=Piecewise[{{{Cos[a+1.2-.48],Sin[a+1.2-.48]},a>0.48}},{Cos[1.2],Sin[1.2]}];
zz[a_]=Piecewise[{{{Cos[a+1.2],Sin[a+1.2]},a<0.72}},{Cos[1.92],Sin[1.92]}];
Manipulate[GraphicsGrid[{{ParametricPlot[{Evaluate[zz[a]],{a^2/2,1/2+Sin[a]}-a/4},{a,0,t},LabelStyle->{Black,Bold,{FontSize->18}},PlotStyle->{{Red,Thick},{Blue,Thick}},Axes->False,Frame->True, ImageSize->{400,300},PlotRange->{{-1,1},{0.2,1.6}},FrameTicks->None,FrameLabel->{x,y},PlotLabel->"The red and blue solutions are never\n in the same position at the same time"]},{ParametricPlot[{Evaluate[z[a]],{a^2/2,1/2+Sin[a]}-a/4},{a,0,t},LabelStyle->{Black,Bold,{FontSize->18}},PlotStyle->{{Red,Thick},{Blue,Thick}},Axes->False,Frame->True, ImageSize->{400,280},PlotRange->{{-1,1},{0.2,1.6}},FrameTicks->None,FrameLabel->{x,y},PlotLabel->"The red solution is now time-shifted"]}}],{t,0.001,1.2,.05}]

Chapter 3

3.3 Examples with Two Pairs of Straight-Line Solutions

A={{1,2},{3,3}};Ei=Eigenvectors[N[A]];Line[{-100Ei[[1]],100Ei[[1]]}];Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"},VectorScale->{Tiny,Small,None},Epilog->{Thick,Red,Line[{-100Ei[[1]],100Ei[[1]]}],Line[{-100Ei[[2]],100Ei[[2]]}]}]

A={{0,-2},{1,-3}};Ei=Eigenvectors[N[A]];Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"},VectorScale->{Tiny,Small,None},Epilog->{Thick,Red,Line[{-100Ei[[1]],100Ei[[1]]}],Line[{-100Ei[[2]],100Ei[[2]]}]}]

A={{0,-1},{2,3}};Ei=Eigenvectors[N[A]];Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"},VectorScale->{Tiny,Small,None},Epilog->{Thick,Red,Line[{-100Ei[[1]],100Ei[[1]]}],Line[{-100Ei[[2]],100Ei[[2]]}]}]

3.4 Examples with Complex Eigenvalues

A={{1,-2},{1,1}};Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];p1=VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,PlotRange->{{-2,2},{-2,2}},VectorScale->{Tiny,Small,None},LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"}];
p2=ParametricPlot[Evaluate[{x[t],y[t]}/.NDSolve[{x'[t]==x[t]-2y[t],y'[t]==x[t]+y[t],x[0]==0,y[0]==.2},{x,y},{t,-10,3}]],{t,-10,3},PlotStyle->{Blue,Thick}];
Show[p1,p2]
Clear[p1,p2];

A={{-1,2},{-2,-2}};Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];p1=VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,PlotRange->{{-2,2},{-2,2}},VectorScale->{Tiny,Small,None},LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"}];
p2=ParametricPlot[Evaluate[{x[t],y[t]}/.NDSolve[{x'[t]==-x[t]+2y[t],y'[t]==-2x[t]-2y[t],x[0]==1.5,y[0]==2},{x,y},{t,-1,10}]],{t,-1,10},PlotStyle->{Blue,Thick}];
Show[p1,p2]

A={{0,2},{-2,0}};Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];p1=VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,PlotRange->{{-2,2},{-2,2}},VectorScale->{Tiny,Small,None},LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"}];
p2=ParametricPlot[Evaluate[{x[t],y[t]}/.NDSolve[{x'[t]==2y[t],y'[t]==-2x[t],x[0]==1,y[0]==1},{x,y},{t,-1,10}]],{t,-1,10},PlotStyle->{Blue,Thick}];
Show[p1,p2]

3.5.1 An Example with a Repeated Eigenvalue: One Pair of Straight-Line Solutions

A={{-5,1},{-1,-3}};Ei=Eigenvectors[N[A]];Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,LabelStyle->{Black,Bold,{FontSize->18}},VectorScale->{Tiny,Small,None},AxesLabel->{"x","y"},Epilog->{Thick,Red,Line[{-100Ei[[1]],100Ei[[1]]}]}]

3.5.2 An Example with a Repeated Eigenvalue: All Solutions Lie on Straight Lines

A={{-2,0},{0,-2}};Ei=Eigenvectors[N[A]];Ev=Eigenvalues[N[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,VectorScale->{Tiny,Small,None},LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"}]

3.5.3 A System with Determinant 0

A={{-2,2},{2,-2}};Ei=N[Eigenvectors[A]];Ev=N[Eigenvalues[A]];Print["The eigenvalues of ",A," are ",Ev[[1]]," and ",Ev[[2]],"."];VectorPlot[If[A.{x,y}=={0,0},{0,0},A.{x,y}/Norm[A.{x,y}]],{x,-2,2},{y,-2,2},Axes->True,VectorScale->{Tiny,Small,None},LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"},Epilog->{Thick,Red,Line[{-100Ei[[1]],100Ei[[1]]}]}]

3.6 A Spring with Variable Damping

Animate[Plot[Evaluate[y[t]/.DSolve[{y''[t]+p y'[t]+4 y[t]==0,y[0]==2,y'[0]==-1},y,t]],{t,0,6},LabelStyle->{Black,Bold,{FontSize->18}},PlotRange->{-2.5,2.1},AxesLabel->{"t","y"},ImageSize->600,PlotLabel->StringForm["y''+ `` y'+4.00 y=0",PaddedForm[p//N,{3,2}]],PlotStyle->{Thick,If[p<4,Red,If[p==4,Green,Blue]]}],{p,0,5,2/100},DisplayAllSteps ->True,AnimationDirection->Forward,AnimationRepetitions->1]

Chapter 5

5.1.1 A Damped Pendulum

A comparison of the direction fields for 	
	y'' (t) + .3 y' (t) + sin (y (t)) = 0 		and 		y'' (t) + .3 y' (t) + y (t) = 0 	
near the origin.

p1=VectorPlot[If[{v,y}=={0,0},{0,0},{v,-y-.3v}/(Norm[{v,-y}])],{y,-3,3},{v,-3,3},VectorColorFunction->Function[x,{Opacity[.75],RGBColor[1,0,0]}],VectorScale->{Tiny,Small,None},Axes->True,AxesLabel->{y,v},PlotLabel->"Damped Pendulum (blue) and its Linearization (red)", LabelStyle->{Black,Bold,{FontSize->18}},PlotRangePadding->.1];
p2=VectorPlot[If[{v,y}=={0,0},{0,0},{v,-Sin[y]-.3v}/(Norm[{v,-Sin[y]}])],{y,-3,3},{v,-3,3},VectorScale->{Tiny,Small,None},PlotRangePadding->.1,Axes->True,VectorColorFunction->Function[x,{Opacity[.8],RGBColor[0,0,1]}]];
Show[p1,p2,ImageSize->600]

A comparison of the phase plots for 
y'' (t) + .3 y' (t) + sin (y (t)) = 0  and  y'' (t) + .3 y' (t) + y (t) = 0  with initial conditions y (0) = 0  and  y' (0) = v0.

pend[v0_]:=NDSolve[{y'[t]==v[t],v'[t]==-Sin[y[t]]-0.3v[t],y[0]==0,v[0]==v0},{y,v},{t,0,50}]
lin[v0_]:=NDSolve[{y'[t]==v[t],v'[t]==-y[t]-0.3v[t],y[0]==0,v[0]==v0},{y,v},{t,0,50}]
pendp[v0_]:=ParametricPlot[Evaluate[{y[t],v[t]}/.pend[v0]],{t,0,50},PlotRange->{{-3,8},{-2,3}},AxesLabel->{y,v},ImageSize->600,Ticks->{Table[k π/2,{k,-2,8}],Automatic},GridLines->{Table[k π/2,{k,-2,8}],Automatic},LabelStyle->{Black,Bold,{FontSize->18}},PlotStyle->Thick]
linp[v0_]:=ParametricPlot[Evaluate[{y[t],v[t]}/.lin[v0]],{t,0,50},PlotRange->{{-3,8},{-2,3}},AxesLabel->{y,v},ImageSize->600,Ticks->{Table[k π/2,{k,-2,8}],Automatic},GridLines->{Table[k π/2,{k,-2,8}],Automatic},LabelStyle->{Black,Bold,{FontSize->18}},PlotStyle->{Opacity[.75],Red}]

Animate[Show[pendp[v0],linp[v0]],{v0,1,2.75,.05},DisplayAllSteps ->True,AnimationDirection->Forward,AnimationRepetitions->1]

5.1.2 Linearization

The module below computes the linearized system and its eigenvalues at a chosen point.

Inputs needed are the right-hand sides of the two differential equations, the two variable names, and the point at which to compute the linearized system.  

A typical example is Linearization[v, -.3 v - Sin[y], y, v, {0, 0}]

Linearization[v, -.3 v - Sin[y], y, v, {π, 0}]

5.1.3 Limitations of Linearization

For the three systems below, the origin is an equilibrium point whose linearization has the eigenvalues ± i .

s0=NDSolve[{x'[t]==y[t],y'[t]==-x[t],x[0]==.1,y[0]==.1},{x,y},{t,0,20}];
ParametricPlot[Evaluate[{x[t],y[t]}/.s0],{t,0,20},PlotRange->All,LabelStyle->{Black,Bold,{FontSize->18}},ImageSize->600,PlotLabel->TraditionalForm[ColumnForm[{x'[t]==y[t],y'[t]==-x[t]}]],Epilog->{Red,AbsolutePointSize[10],Point[{.1,.1}]}]

s1=NDSolve[{x'[t]==y[t]-(x[t]^2+y[t]^2)x[t],y'[t]==-x[t]-(x[t]^2-y[t]^2)y[t],x[0]==.1,y[0]==.1},{x,y},{t,0,50}];
ParametricPlot[Evaluate[{x[t],y[t]}/.s1],{t,0,50},PlotRange->All,LabelStyle->{Black,Bold,{FontSize->18}},ImageSize->600,PlotLabel->TraditionalForm[ColumnForm[{x'[t]==y[t]-(x[t]^2+y[t]^2)x[t],y'[t]==-x[t]-(x[t]^2-y[t]^2)y[t]}]],Epilog->{Red,AbsolutePointSize[10],Point[{.1,.1}]}]

s2=NDSolve[{x'[t]==y[t]+(x[t]^2+y[t]^2)x[t],y'[t]==-x[t]+(x[t]^2-y[t]^2)y[t],x[0]==.1,y[0]==.1},{x,y},{t,0,95}];
ParametricPlot[Evaluate[{x[t],y[t]}/.s2],{t,0,95},PlotRange->All,LabelStyle->{Black,Bold,{FontSize->18}},ImageSize->600,PlotLabel->TraditionalForm[ColumnForm[{x'[t]==y[t]+(x[t]^2+y[t]^2)x[t],y'[t]==-x[t]+(x[t]^2-y[t]^2)y[t]}]],Epilog->{Red,AbsolutePointSize[10],Point[{.1,.1}]}]

5.1.4 Separatrices - A Toy Example

The system
		x' = -x
		 y' = -4 x^3 + y
is non-linear, but partially decoupled; therefore we can find its solutions explicitly.
The only equilibrium point is (0,0); its linearization yields a saddle:

Linearization[-x, -4x^3+y, x, y, {0, 0}]

Below is the direction field of the system; the two stable separatrices are displayed in red, the two unstable separatrices in blue. 
Notice the crtitical role played by the stable separatrices in determining the long-term behavior of the solutions.

p1=VectorPlot[If[{-x,-4x^3+y}=={0,0},{0,0},{-x,-4x^3+y}/Norm[{-x,-4x^3+y}]],{x,-1.5,1.5},{y,-2,2},Axes->True,LabelStyle->{Black,Bold,{FontSize->18}},AxesLabel->{"x","y"},VectorScale->{Tiny,Small,None},Epilog->{AbsoluteThickness[4],Blue,Line[{{0,-2},{0,2}}]},AspectRatio->1,ImageSize->600,PlotRangePadding->0.125];
p2=Plot[x^3,{x,-1.5,1.5},PlotStyle->{AbsoluteThickness[4],Red},PlotRange->{-2,2}];
Show[p1,p2]

5.1.5 Competing Species - Basins of Attraction

The system 
		 x'= 2x (1 - x/2) - xy
			 y'= 3y (1 - y/3) - 2xy
describes  the interaction between two competing species.
The model has four equilibrium points:

Solve[{0==2.x[t](1.-x[t]/2.)-x[t] y[t],0==3.y[t](1.-y[t]/3.)-2.x[t] y[t]},{x[t],y[t]}]

Linearization classifies the equilibrium points as follows :
	 	 (0,0)  is a source - corresponding to the extinction of both species.
		 (2,0) and (0,3) are sinks - corresponding to the extinction of one of the two species.
		 (1,1) is a saddle - corresponding to co-existence of the two species.
	
Below we color initial conditions in the (x,y)-plane according to the fate of solutions "in the long run". Note that all solutions in the first quadrant approach one of the four equilibrium points. 
If a solution approaches 	(0,3) 	the initial condition is colored		blue
				(2,0)						red
				(0,0)						black
				(1,1)						white
The dividing line between the blue and the red region is (an approximation of) the phase portrait of the two stable separatrices.

The last picture takes several hours of computing time.

sol[a_,b_]:=NDSolve[{x'[t]==2.x[t](1.-x[t]/2.)-x[t] y[t],y'[t]==3.y[t](1.-y[t]/3.)-2.x[t] y[t],x[0]==a,y[0]==b},{x,y},{t,0,5000}]
Do[tu=TimeUsed[];rlist=Table[If[Flatten[Chop[{x[t],y[t]}/.sol[a,b]/.t->5000],1]=={0,3.},{0,0,1},If[Flatten[Chop[{x[t],y[t]}/.sol[a,b]/.t->5000],1]=={2.,0},{1,0,0},
If[Flatten[Chop[{x[t],y[t]}/.sol[a,b]/.t->5000],1]=={1.,1.},{1,1,1},{0,0,0}]]],{b,3,0,-10.^-k},{a,0,3,10.^-k}];Print[ArrayPlot[1-rlist,DataRange->{{0,3},{0,3}},Frame->True,FrameLabel->{"y","x"}, PlotLabel->StringForm["step size = `` \t(`` initial conditions)",10.^(-k),(3 10^k+1)^2],FrameTicks->{{Automatic,None},{Automatic,None}},ImageSize->700,LabelStyle->{Black,Bold,{FontSize->18}}]];Print[Ceiling[tu=TimeUsed[]-tu]," sec."],{k,0,3}]

5.2.1 Nullcline Analysis

The command below produces a nullcline plot for a system of differential equations. 
Inputs needed are the right-hand sides of the two differential equations, and the ranges for the two variables. 
Here is a typical example:
 	NullClinePlot[.25 x^2 + .5 y^2 - 3, (x - .75)^3 + 2 y^2 - 5, {x, -4,  4}, {y, -4, 4}]

NullClinePlot[.25x^2+.5y^2-3,(x-.75)^3+2y^2-5,{x,-4,4},{y,-4,4}]

5.2.2 Competing Species - Linearization and Nullcline Analysis

The system 
	x' = x (100 - x - 2 y)
	y' = y (150 - x - 6 y)
describes two competing species. The system has 4 equilibrium points.

The equilibrium point corresponding to coexistence of the two species is a sink.
The two equilibrium points corresponding to extinction of one of the two species are saddles;  their stable separatrices lie on the coordinate axes.

Linearization[x(100-x-2y),y(150-x-6y), x, y, {0,0}]

Linearization[x(100-x-2y),y(150-x-6y), x, y, {100,0}]

Linearization[x(100-x-2y),y(150-x-6y), x, y, {0,25}]

Solve[{100-x-2y==0,150-x-6y==0}]
Linearization[x(100-x-2y),y(150-x-6y), x, y, {75,25/2}]

NullClinePlot[x(100-x-2y),y(150-x-6y),{x,-20,200},{y,-5,55}]

Chapter 6

6.2 The Unit Step Function

Plot[UnitStep[t-3],{t,0,7},PlotStyle->{AbsoluteThickness[3],Red},ImageSize->600]

We can use unit step functions to represent piecewise defined functions:

Clear[t,s,sol];
f[t_]=t^2+UnitStep[t-1](t-t^2)+UnitStep[t-2]((3-t/2)-t)+UnitStep[t-4](-(3-t/2));
Plot[f[t],{t,0,6},PlotStyle->{AbsoluteThickness[3],Red},ImageSize->600]
PiecewiseExpand[f[t]]

Let' s solve a differential equation with this function as the right - hand side. 

Note that Mathematica displays the unit step function as θ(t).

eq=y'[t]+3y[t]==f[t]
s1=LaplaceTransform[eq,t,s]

s2=Solve[s1,LaplaceTransform[y[t],t,s]]/.y[0]->1/2;
s3=s2[[1,1,2]]

sol[t_]=InverseLaplaceTransform[s3,s,t]

Here is a plot of the solution and its derivative:

GraphicsGrid[{{Plot[sol[t],{t,0,5},PlotStyle->{Thick, Blue},AxesLabel->{t,y},BaseStyle->{14,Bold},PlotRange->All,PlotLabel->"y(t)"]},{
Plot[sol'[t],{t,0,5},PlotStyle->{Thick, Orange},AxesLabel->{t,y},BaseStyle->{14,Bold},PlotRange->All,PlotLabel->"y'(t)"]}},ImageSize->400]

6.3 Interference and Resonance

We consider an undamped second - order system subject to a sinosoidal input force:

	y''[t]+((Subscript[ω, 0])^2)y[t]=A Sin[(Subscript[ω,  ])t]

If the period ω of the input function is close to the natural period (Subscript[ω, 0]) of the system (in our scenario below this natural period is (Subscript[ω, 0])=4), the solution shows the phenomenon called interference; if the input period and the natural periods coincide, we obtain resonance, i.e the solution becomes unbounded.

These phenomena only occur if the system is undamped.

ω0=4;ω=4.0005;
Clear[s,t,sol];
lt=LaplaceTransform[y''[t]+ω0^2y[t]==125Sin[ω t],t,s]

lt2=Solve[lt,LaplaceTransform[y[t],t,s]]/.{y[0]->1000,y'[0]->3000}

sol[t_]=Chop[FullSimplify[ComplexExpand[InverseLaplaceTransform[lt2[[1,1,2]],s,t]]]]

We can now "listen" to the solution:

Play[sol[500t],{t,0,2},SampleRate->10000]