Ordinary Differential Equations
and
Dynamical Systems
with Mathematica

©2001-2020 by Gerald Teschl  <mailto:Gerald.Teschl@univie.ac.at>, <http://www.mat.univie.ac.at/~gerald/>

First order equations

The command VectorPlot will illustrate the vector field corresponding to a given equation.

VectorPlot[{1,x Sin[t]},{t,0,2 π},{x,0,6}]

An alternative is the StreamPlot command

StreamPlot[{1,(x^2)Sin[t]},{t,0,2 π},{x,0,6}]

Now here is the general solution of this equation.

DSolve[x'[t]==x[t] Sin[t],x[t],t]

And here is the solution of the initial value problem.

DSolve[{x'[t]==x[t] Sin[t],x[0]==1},x[t],t]

Here is how it looks like.

Plot[x[t]/.%,{t,0,2 π}]

Only one solution is found here:

DSolve[{x'[t]==(Sqrt[x[t] ]),x[0]==0},x[t],t]

The first "solution" is no solution!

DSolve[{x'[t]==(Sqrt[x[t] ]),x[0]==x0},x[t],t]//Simplify

It satisfies x'[t]==-(Sqrt[x[t] ])

Qualitative analysis

Riccati equation

NDSolve[{x'[t]==x[t]^2-t^2,x[0]==1},x[t],{t,-2,2}]

An overview can be obtained using StreamPlot:

StreamPlot[{1,(x^2)-(t^2)},{t,-4,4},{x,-3,3}]

Alternatively one can use the PhasePlot command from section "Planar dynamical systems" below.

PhasePlot[{1,y^2-x^2},{{0,-2},{1,-1.5},{0,-1.25},{0,-0.75},{0,-0.676},{0,-0.5},{0,-0.25},{0,0},{0,0.1},{0,0.25},{0,0.5},{0,0.676},{0,0.75},{-1,1.5},{0,2},{1,1},{-1,-1},{1.5,1.5},{-1.5,-1.5},{2,2},{-2,-2},{-1.8,2.3},{1.8,-2.3},{3,2.1},{-3,-2.1}},4,PlotRange->{{-4,4},{-3,3}},Axes->True,AxesLabel->{"t","x"}]

xp[t_]=t; ym[t_]=(Sqrt[2+t^2]);
x1[t_]= x[t] /.Flatten[NDSolve[{x'[t]==x[t]^2-t^2,x[0.5]==xp[0.5]},x[t],{t,0,1.2}]];
x2[t_]= x[t] /.Flatten[NDSolve[{x'[t]==x[t]^2-t^2,x[0.5]==ym[0.5]},x[t],{t,0,1.2}]];

gr1=Plot[{xp[t],ym[t]},{t,0,1.2`},Filling->{1->{{2},GrayLevel[0.8`]}}];
gr2=Plot[{x1[t],x2[t]},{t,0,1.2`}];Show[gr1,gr2,Ticks->{{{0.5,"T=0.5"},1,1.5},{0,{x1[0],"a(T)"},{x2[0],"b(T)"},1,2}},PlotRange->{0,2.3}]

Logistic equation with periodic harvesting

x[t] /. DSolve[{x'[t]== x[t](1-x[t]),x[0]==x0},x[t],t][[1]] /. t->1

p[x0_,h_]:= Block[{sol,x,t},sol=x /.NDSolve[{x'[t]==x[t](1-x[t])-h(1+Sin[2π t]),x[0]==x0},x,{t,0,1}][[1]];
sol[1] 
]

GraphicsRow[{Plot[{p[x,0],x},{x,0,1},Ticks->{{0,1/2,1},{0,1/2,1}},PlotRange->{-0.7,1},PlotLabel->"h=0",AspectRatio->1],Plot[{p[x,0.1],x},{x,0,1},Ticks->{{0,1/2,1},{0,1/2,1}},PlotRange->{-0.7,1},PlotLabel->"h=0.1",AspectRatio->1],Plot[{p[x,0.249],x},{x,0,1},Ticks->{{0,1/2,1},{0,1/2,1}},PlotRange->{-0.7,1},PlotLabel->"h=0.249",AspectRatio->1],Plot[{p[x,0.5],x},{x,0,1},Ticks->{{0,1/2,1},{0,1/2,1}},PlotRange->{-0.7,1},PlotLabel->"h=0.5",AspectRatio->1]}]

Export["poincaremap.jpg",%]

PhasePlot[{1,y(1-y)-0.2(1+Sin[2π x])},{{0,0.2},{0,0.3},{0,0.6},{0,1},{0,2},{0,6}},4,PlotRange->{{0,4},{0,2}},Ticks->False]

StreamPlot[{1,x(1-x)-0.2(1+Sin[2π t])},{t,0,2 π},{x,-1,2}]

SIS model

sol=DSolve[{i'[t]== β ((i[t])/N0)(N0-i[t]) - γ i[t],i[0]==i0},i[t],t][[1]]//FullSimplify

i[t] /. sol /.N0->N1/(1-(γ/β))//FullSimplify

Limit[i[t] /. sol , γ->β]

Manipulate[Plot[{(( i0 (1-γ))/( i0+(E^(-t(1- γ))) ((1-γ)-i0))),(1-γ)},{t,0,5},PlotRange->{-0.1,1}],{i0,0.1,1},{γ,0,3}]

The contraction principle

Iteration of a function:

ShowWeb[K_,xstart_,nmax_]:=Block[{x,xmin,xmax,graph,web},
x[0]:= xstart;
x[n_]:=x[n]=K[x[n-1]];
If[x[0]==x[1],Return[Print[Fixed Point!]]];
web=Flatten[Table[{{x[n],x[n]},{x[n],x[n+1]}},{n,0,nmax}],1];
xmax=Max[web]; xmin=Min[web];
graph= Plot[{K[x],x},{x,xmin,xmax},DisplayFunction->Identity];
Show[graph,Graphics[Line[web]],DisplayFunction->$DisplayFunction]
] /; (IntegerQ[nmax] && Positive[nmax])

ShowWeb[(1-((#)^2)/2)&,0,12]

ShowWeb[-2ArcTan[#]&,1,12]

Animation

Manipulate[ShowWeb[(1-((#)^2)/2)&,0,nmax],{nmax,1,12,1}]

Picard iteration

f[t_,x_]=((1 )/(1-x));x0=0;

x[0,t_]=x0;
x[n_,t_]:=x0+(∫_0^t)f[s,x[n-1,s]]\[DifferentialD]s;

Series[x[2,t],{t,0,2}]

Series[x[t] /. DSolve[{x'[t]==f[t,x[t]],x[0]==x0},x[t],t][[1]],{t,0,2}]//FullSimplify

Numerics

Clear[f,sol,t0,t1,nmax]

f[t_,x_]=x;t0=0;t1=10;x0=1;nmax=20;

sol[t_]=x[t]/.DSolve[{x'[t]==f[t,x[t]],x[t0]==x0},x[t],t][[1]]

Euler

EulerSolve[f_,t0_,t1_,x0_,nmax_]:=Block[{tmp={{t0,N[x0]}},h=((t1-t0)/nmax),tt,xx},Do[{tt,xx}=tmp[[-1]];xx=xx+ f[tt,xx]h;AppendTo[tmp,{tt+h,xx}],{n,1,nmax}];
tmp];

sole=EulerSolve[f,t0,t1,x0,nmax];

N[{sole[[-1,2]],sol[t1]}]

Show[ListPlot[sole,PlotStyle->{ColorData[1,2],PointSize[0.012]}],Plot[sol[x],{x,t0,t1}]]

Heun

HeunSolve[f_,t0_,t1_,x0_,nmax_]:=Block[{tmp={{t0,N[x0]}},h=((t1-t0)/nmax),tt,xx,k1},Do[{tt,xx}=tmp[[-1]];k1=f[tt,xx];xx=xx+( k1+f[tt+h,xx+k1 h])(h/2);AppendTo[tmp,{tt+h,xx}],{n,1,nmax}];
tmp];

solh=HeunSolve[f,t0,t1,x0,nmax];

N[{solh[[-1,2]],sol[t1]}]

Show[ListPlot[solh,PlotStyle->{ColorData[1,3],PointSize[0.012]}],Plot[sol[x],{x,t0,t1}]]

Runge-Kutta

RK4Solve[f_,t0_,t1_,x0_,nmax_]:=Block[{tmp={{t0,N[x0]}},h=((t1-t0)/nmax),tt,xx,k1,k2,k3,k4},Do[{tt,xx}=tmp[[-1]];k1=f[tt,xx];
k2=f[tt+(h/2),xx+(h/2)k1];
k3=f[tt+(h/2),xx+(h/2)k2];
k4=f[tt+h,xx+h k3];
xx=xx+ (h/6)(k1+2k2+2k3+k4);AppendTo[tmp,{tt+h,xx}],{n,1,nmax}];
tmp];

solrk=RK4Solve[f,t0,t1,x0,nmax];

N[{solrk[[-1,2]],sol[t1]}]

Show[ListPlot[solrk,PlotStyle->{ColorData[1,4],PointSize[0.012]}],Plot[sol[x],{x,t0,t1}]]

Error comparison (equal step size)

sole=EulerSolve[f,t0,t1,x0,nmax];solh=HeunSolve[f,t0,t1,x0,nmax];solrk=RK4Solve[f,t0,t1,x0,nmax];

N[sol[t1]-{sole[[-1,2]],solh[[-1,2]],solrk[[-1,2]]},2]

Error comparison (equal number of evaluations of f)

sole=EulerSolve[f,t0,t1,x0,4nmax];solh=HeunSolve[f,t0,t1,x0,2nmax];solrk=RK4Solve[f,t0,t1,x0,nmax];

N[sol[t1]-{sole[[-1,2]],solh[[-1,2]],solrk[[-1,2]]},2]

Visualize the data

Clear[nmax];

f[t_,x_]=(-1/10+Sin[t])x;t0=0;t1=10;x0=1;

f[t_,x_]=-t (x^2);t0=0;t1=10;x0=1;

f[t_,x_]=(t^2)- x;t0=0;t1=20;x0=1;

sol[t_]=x[t]/.DSolve[{x'[t]==f[t,x[t]],x[t0]==x0},x[t],t][[1]]

MyOpacity[x_]=If[x,Opacity[1],Opacity[0]];Attributes[MyOpacity]={Listable};

Manipulate[Show[Plot[sol[t],{t,t0,t1},PlotRange->All,PlotStyle->GrayLevel[0.5]],ListLinePlot[{EulerSolve[f,t0,t1,x0,nmax],HeunSolve[f,t0,t1,x0,nmax],RK4Solve[f,t0,t1,x0,nmax]},PlotRange->All,PlotMarkers->Automatic,PlotStyle->MyOpacity[{e,h,r}]]],{{nmax,10},1,50,1},Row[{Control[{{e,True,"Euler"},{True,False}}],Control[{{h,False,"Heun"},{True,False}}],Control[{{r,True,"RK4"},{True,False}}]},Spacer[20]]]

Linear systems

The Jordan decomposition of a matrix

A=(({{-11, -35, -24}, {-1, -1, -2}, {8, 22, 17}}));

{U,J}=JordanDecomposition[A];
J // MatrixForm

A==U.J.Inverse[U]

MatrixExp[t J] //MatrixForm

The real Joran decomposition of a matrix

A=(({{2, -1, -1, 0}, {1, 0, 0, -1}, {1, 0, 2, -1}, {0, 1, 1, 0}}));

{U,J}=JordanDecomposition[A];
J // MatrixForm

U // MatrixForm

TU=Transpose[U];
V=Transpose[{Re[TU[[1]]],Im[TU[[1]]],Re[TU[[2]]],Im[TU[[2]]]}];
V// MatrixForm

Jr=Inverse[V].A.V;
Jr//MatrixForm

MatrixExp[t Jr]//FullSimplify//MatrixForm

Clear[U,V,TU,Jr]

The logarithm of a Jordan block

Nil[n_Integer,m_Integer]:=Table[If[j-i==m,1,0],{i,1,n},{j,1,n}];
Ln[n_Integer]:=Log[α]Nil[n,0]-(∑_(j=1)^n)((((-1))^j)/(j (α^j)))Nil[n,j];

Ln[4] //MatrixForm

MatrixExp[Ln[4]] //Simplify //MatrixForm

Nil[n_]:=Table[If[j-i==1,1,0],{i,n},{j,n}];
Nil[5]//MatrixForm

U[n_]:=Table[((( j-1)!StirlingS1[k-1,j-1])/((k-1)!)),{j,n},{k,n}];
V[n_]:=Table[((( j-1)!StirlingS2[k-1,j-1])/((k-1)!)),{j,n},{k,n}];

V[5]==Inverse[U[5]]

V[5].MatrixExp[Nil[5]].U[5]==IdentityMatrix[5]+Nil[5]

The Routh-Hurwitz criterion

RouthHurwitzQ[p_,x_]:=Block[{cl=Reverse[CoefficientList[p,x]],l},
l=Length[cl];Reduce[Table[
Det[Table[If[0≤2i-j<l,cl[[2i-j+1]],0],{i,1,k},{j,1,k}]]>0,
{k,l-1}]]
] /; PolynomialQ[p,x];
RouthHurwitzQ[m_?MatrixQ]:=RouthHurwitzQ[Det[IdentityMatrix[Dimensions[m][[1]]]x-m],x];

RouthHurwitzQ[(({{-1, 1}, {1, -2}}))]

The following function computes the first column of the Routh array. The number of sing flips gives the number of eigenvalues with positive real part.

RouthArray[p_,x_]:=Block[{a=Reverse[CoefficientList[p,x]],n,m,RA},
n=Length[a]-1; 
m=Ceiling[((n+1)/2)];
RA=Table[a[[2j-1]],{j,m}];
RA={RA,Table[If[2j≤n+1,a[[2j]],0],{j,m}]};
Do[
RA=Append[RA,Table[If[j<m,((RA[[l-1,1]]RA[[l-2,j+1]]-RA[[l-2,1]]RA[[l-1,j+1]])/(RA[[l-1,1]])),0]
,{j,m}]],
{l,3,n+1}];
Simplify[Transpose[RA][[1]]]
]/; PolynomialQ[p,x];
RouthArray[m_?MatrixQ]:=RouthArray[Det[IdentityMatrix[Dimensions[m][[1]]]x-m],x];

RouthHurwitzQ[(({{-1, 1}, {1, 0}}))]

RouthArray[(({{-1, 1}, {1, 0}}))]

Eigenvalues[(({{-1, 1}, {1, 0}}))]

Linear systems - warning

Mathematica up to version 9.0 can give wrong solutions for linear systems!

A=(({{-24/10, 73/100}, {49/100, -58/100}}));x0={0,0};
eq=D[{x[t],y[t]},t]==A . {x[t],y[t]}+{0,f};

The solution via the variation of constants formula

sol1=MatrixExp[A t].x0+(∫_0^t)MatrixExp[A(t-s)].{0,f}\[DifferentialD]s//N//Expand//Simplify

The solution via DSolve

sol2={x[t],y[t]} /.DSolve[{ eq,  {x[0], y[0]}== x0},{x[t],y[t]},t][[1]]//N//Expand//Simplify

So far so good

sol1-sol2//Simplify//Chop

Specifying the coefficient matrix A as decimal values gives a wrong solution

sol3={x[t],y[t]} /.DSolve[{ N[eq],   {x[0], y[0]}== x0},{x[t],y[t]},t][[1]]//N//Expand//Simplify

sol1-sol3//Simplify//Chop

RLC circuit

sol=Simplify[DSolve[{x''[t]+4x'[t] +x[t]==0,x[0]==1,x'[0]==1},x[t],t][[1]]]

Plot[x[t]/.sol,{t,0,12},Ticks->None,PlotLabel->"over damping"]

sol=Simplify[DSolve[{x''[t]+2x'[t] +x[t]==0,x[0]==1,x'[0]==1},x[t],t][[1]]]

Plot[x[t]/.sol,{t,0,12},Ticks->None,PlotLabel->"critical damping"]

sol=Simplify[DSolve[{x''[t]+x'[t] /2+x[t]==0,x[0]==1,x'[0]==1},x[t],t][[1]]]

Plot[x[t]/.sol,{t,0,12},Ticks->None,PlotLabel->"under damping"]

Mathieu

MathieuEq[δ_,ϵ_]={x'[t]==v[t],v'[t]==- (δ^2)(1+ϵ Cos[t])x[t]};

spur[δ_?NumericQ,ϵ_?NumericQ]:=Block[{t,c,sp, ivp},
ivp=Flatten[{MathieuEq[δ,ϵ],x[0]==1,v[0]==0}];
c=(x[t] /.Flatten[NDSolve[ivp,{x[t],v[t]},{t,0,2π}]]) /. t->2π;
ivp=Flatten[{MathieuEq[δ,ϵ],x[0]==0,v[0]==1}];
sp=(v[t] /.Flatten[NDSolve[ivp,{x[t],v[t]},{t,0,2π}]]) /. t->2π;
((c+sp)/2)
];

ContourPlot[((spur[ω,ϵ])^2),{ω,0,3},{ϵ,-1.5`,1.5`},PlotPoints->70,AspectRatio->1,ContourShading->{None,Gray},Contours->{0.999},ClippingStyle->Automatic,ContourStyle->Black,FrameLabel->Automatic]

Export["ode3_7.eps",%];

Complex systems, power series

This section requires the package DiscreteMath'DiffEqs' available from
	 http://www.mat.univie.ac.at/~gerald/ftp/book-jac/DiffEqs.m
It needs to be  stored in the Applications/DiscreteMath subfolder of your Mathematica folder. On multi-user systems, you can install an add-on either in the central Mathematica directory (provided you have access to it) or else in your individual user's Mathematica directory (usually ~/.Mathematica/). (You can also put them  somewhere else, as long as you change the $Path variable in Mathematica such that the files can be found by Mathematica or evaluate the packages before using this notebook.)

FuchsParameters[p_,q_,z_,∞]:=Block[{pp0,qq0,α,α1,α2},
pp0=2-Limit[z p,z->∞];
qq0=Limit[(z^2)q,z->∞];
{α1,α2}= α /.Solve[(α^2)+(pp0-1)α+qq0==0,α];
{pp0,qq0,α1,α2}];
FuchsParameters[p_,q_,z_,z0_]:=Block[{pp0,qq0,α,α1,α2},
pp0=Limit[(z-z0)p,z->z0];
qq0=Limit[(((z-z0))^2)q,z->z0];
{α1,α2}= α /.Solve[(α^2)+(pp0-1)α+qq0==0,α];
{pp0,qq0,α1,α2}];

Linear difference equation

Needs["DiscreteMath`DiffEqs`"]

xh[n_]:=x0 SymbProduct[f[j] ,{j,0,n-1}];
x[n_]:=xh[n]+xh[n]SymbSum[((g[k])/(xh[k+1])) ,{k,0,n-1}];

x[n+1]-f[n]x[n]==g[n]//SplitAll

Riccati

Needs["DiscreteMath`DiffEqs`"]

RiccatiEqn[f_]:=(∂_z) f  + (f^2)-z;

f[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[0]:=h0;h[j_]:=0 /; (j<0);

eqn=RiccatiEqn[f[z]]

FactorSum[SplitSum[ShiftPowers[ExpandSum[CauchyProduct[eqn]],z],0]] //Simplify

Linear

Needs["DiscreteMath`DiffEqs`"]

MySeries[f_,{z_,z0_}]:= Block[{ff},
ff[j_]=FullSimplify[SeriesCoefficient[f,{z,z0,j}],j∈Integers&&j≥0];
SymbSum[ff[j] (z^j),{j,0,∞}]
];

g[z_]=MySeries[(1/(1-z)),{z,0}]

LinearEqn[f_]:=(∂_z) f  + f-g[z];

f[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[0]:=h0;h[j_]:=0 /; (j<0);

eqn=LinearEqn[f[z]]

FactorSum[SplitSum[ShiftPowers[ExpandSum[eqn],z],0]] //Simplify

RSolve[-1+h[j]+(1+j) h[1+j]==0,h[j],j]//FullSimplify[#,j∈Integers&&j≥0]&

Linear 2

Needs["DiscreteMath`DiffEqs`"]

LinearEqn[f_]:=(∂_(z,z)) f  + (1-(z^2))f;

f[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[j_]:=0 /; (j<0);

eqn=LinearEqn[f[z]]

FactorSum[SplitSum[ShiftPowers[ExpandSum[eqn],z],0]] //FullSimplify

Clear[h]

RSolve[{h[-2+j]-h[j]-(1+j) (2+j) h[2+j]==0,h[-2]==0,h[-1]==0,h[0]==1,h[1]==0},h[j],j]

RSolve[{h[-1+j]-h[j]-(1+2j) (2+2j) h[1+j]==0,h[0]==1,h[-1]==0},h[j],j]

Bessel

Needs["DiscreteMath`DiffEqs`"]

FuchsParameters[(1/z),(((z^2)- (ν^2))/(z^2)),z,0]

BesselEqn[f_]:=(z^2)(∂_(z,z)) f+z (∂_z) f  + ((z^2)- (ν^2) )f;

u[z_]:=SymbSum[h[j] (z^(j+ν)),{j,0,∞}];
h[0]:=1; h[j_]:=0 /; (j<0);

eqn=BesselEqn[u[z]]

ShiftPowers[ExpandSum[ (z^(-ν)) eqn],z]

SplitSum[%,0] // FactorSum //Simplify

hh[j_]=hh[j] /. Flatten[RSolve[{hh[j-1]+4j (j+ ν) hh[j]==0,hh[0]==1},hh[j],j]]

v[z_]:=SymbSum[g[j] (z^(j-ν))+ c Log[z] h[j] (z^(j+ν)),{j,0,∞}];
g[0]:=1; g[j_]:=0 /; (j<0);

eqn=BesselEqn[v[z]] /. Log[z]->0

FactorSum[%]//Simplify

Unprotect[IntegerQ];
IntegerQ[2ν]=True;

ShiftPowers[ExpandSum[ (z^ν) eqn],z]

SplitSum[%,0] // FactorSum //Simplify

RSolve[{gg[-1+j]+4 j (j+ν) gg[j]+2 c (2 j+ν) hh[j]==0,gg[1]==1},gg[j],j]

Plot[{BesselJ[0,z],BesselY[0,z]},{z,0,12}]

z(BesselJ[n,z]D[BesselY[n,z],z] -D[BesselJ[n,z],z]BesselY[n,z] )//FullSimplify

Sum[2k+ν,{k,1,n-1}]

Legendre

Needs["DiscreteMath`DiffEqs`"]

FuchsParameters[((-2z)/(1-(z^2))),((n(n+1))/(1-(z^2))),z,1]

FuchsParameters[((-2z)/(1-(z^2))),((n(n+1))/(1-(z^2))),z,-1]

FuchsParameters[((-2z)/(1-(z^2))),((n(n+1))/(1-(z^2))),z,∞]

LegendreEqn[f_]:=(1-(z^2)) (∂_(z,z)) f-2z (∂_z) f  + n(n+1)f;

u[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[0]:=1; h[j_]:=0 /; (j<0);

eqn=LegendreEqn[u[z]]

ShiftPowers[ExpandSum[  eqn],z]

SplitSum[%,0] // FactorSum //FullSimplify

Hermite

Needs["DiscreteMath`DiffEqs`"]

HermiteEqn[f_]:= (∂_(z,z)) f-2z (∂_z) f  + 2 n f;

u[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
 h[j_]:=0 /; (j<0);

eqn=HermiteEqn[u[z]]

ShiftPowers[ExpandSum[  eqn],z]

SplitSum[%,0] // FactorSum //FullSimplify

h[j] /. RSolve[{2 (-j+n) h[j]+(1+j) (2+j) h[2+j]==0,h[0]==1,h[1]==0},h[j],j][[1]]

h[j] /. RSolve[{2 (-j+n) h[j]+(1+j) (2+j) h[2+j]==0,h[0]==0,h[1]==1},h[j],j][[1]]

Hypergeometric

Needs["DiscreteMath`DiffEqs`"]

FuchsParameters[((c-(1+a+b)z)/(z(1-z))),((-a b)/(z(1-z))),z,0]

FuchsParameters[((c-(1+a+b)z)/(z(1-z))),((-a b)/(z(1-z))),z,1]

FuchsParameters[((c-(1+a+b)z)/(z(1-z))),((-a b)/(z(1-z))),z,∞]

HypergeometricEqn[f_]:=z(1-z)(∂_(z,z)) f+(c-(1+a+b)z) (∂_z) f  - a b f;

u[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[0]:=1; h[j_]:=0 /; (j<0);

eqn=HypergeometricEqn[u[z]]

ShiftPowers[ExpandSum[ eqn],z]

SplitSum[%,0] // FactorSum //Simplify

Confluent Hypergeometric

Needs["DiscreteMath`DiffEqs`"]

FuchsParameters[((c-z)/z),((-a )/z),z,0]

ConfluentHypergeometricEqn[f_]:=z(∂_(z,z)) f+(c-z) (∂_z) f  - a f;

u[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];
h[0]:=1; h[j_]:=0 /; (j<0);

eqn=ConfluentHypergeometricEqn[u[z]]

ShiftPowers[ExpandSum[ eqn],z]

SplitSum[%,0] // FactorSum //Simplify

Riemann equation

ζ[z_]=(((z-z0)(z1-z2))/((z-z2)(z1-z0)));

z/.Solve[ζ[z]==zeta,z][[1]]//Together

FuchsParameters[(p0/z)+(p1/(z-1)),(q0/(z^2))+(r0/z)+(q1/(((z-1))^2))-(r0/(z-1)),z,0]//FullSimplify

FuchsParameters[(p0/z)+(p1/(z-1)),(q0/(z^2))+(r0/z)+(q1/(((z-1))^2))-(r0/(z-1)),z,1]//FullSimplify

FuchsParameters[(p0/z)+(p1/(z-1)),(q0/(z^2))+(r0/z)+(q1/(((z-1))^2))-(r0/(z-1)),z,∞]//FullSimplify

Sturm-Liouville

q[x_]=0;

θ[λ_]:=Block[{x,f},f /.NDSolve[{f'[x]==((Cos[f[x]])^2)+(λ-q[x])((Sin[f[x]])^2),f[0]==0},f,{x,0,π}][[1]]]

Plot[θ[λ][π]/π,{λ,-5,10},GridLines->{{},{1,2,3}},AxesLabel->{"λ","(Subscript[θ, a])(λ,b)/π"}]

Export["ode5_1.eps",%];

Plot[Table[θ[λ][x]/π,{λ,-1,5}],{x,0,π},GridLines->{{},{1,2}},AxesLabel->{"x","(Subscript[θ, a])(λ,x)/π"}]

Export["ode5_2.eps",%];

Prüfer transform

u[x_]=((ρ[x]Sin[θ[x]])/(Sqrt[h[x]]));
pup[x_]=(Sqrt[h[x]])ρ[x]Cos[θ[x]];

sol=(FullSimplify[Solve[{D[pup[x],x] == (q[x]-λ r[x])u[x] ,p[x]D[u[x],x]== pup[x]},{θ'[x],ρ'[x]}][[1]] ,h[x]>0] /. ((Cos[θ_])^2)->1-((Sin[θ])^2)//Expand)

Unmodified Prüfer equation

h[x_]=1;

θ'[x] /. sol//Simplify

((ρ'[x] )/(ρ[x]))/. sol//Expand

High energy asymptotics

h[x_]=(Sqrt[λ r[x] p[x]]);

θ'[x] /. sol//PowerExpand//Expand

((ρ'[x] )/(ρ[x]))/. sol//Expand//PowerExpand

Bessel equation

p[x_]=1;r[x_]=0;h[x_]=(Sqrt[-q[x]]);q[x_]=-1-(((1/4)-(ν^2))/(x^2));

θ'[x] /. sol/. θ[x]->θ//Simplify

Series[%,{x,∞,2}]

((ρ'[x] )/(ρ[x]))/. sol/. θ[x]->θ//FullSimplify

Series[%,{x,∞,2}]//Simplify

Periodic Sturm-Liouville (Hille) equation

DSolve[x''[t]+x[t]==Cos[α t],x[t],t]//FullSimplify

DSolve[x''[t]+x[t]==Cos[t],x[t],t]//FullSimplify

Quantum mechanics of an electron in a crystal

Φ[z_,t_]:=(({{Cos[(Sqrt[z])t], ((Sin[(Sqrt[z])t])/(Sqrt[z]))}, {-(Sqrt[z])Sin[(Sqrt[z])t], Cos[(Sqrt[z])t]}}));

Φ[z,1]//MatrixForm

FullSimplify[(1/2)Tr[Φ[z,1]]]

M[z_]=FullSimplify[Φ[z,(1/2)].Φ[z-q0,(1/2)]];

Δ[z_]=FullSimplify[(1/2)Tr[M[z]]]

Δ[z]/. q0->z-(((a-ϵ))^2)/.z->(((a+ϵ))^2)//FullSimplify

q0=50;Plot[{Δ[z],1,-1},{z,0,200},PlotRange->{-1.5,1.5},AxesLabel->{"z","Δ(z)"}]

Export["ode5_3.eps",%];

ρp[z_]=If[Abs[Δ[z]+(Sqrt[((Δ[z])^2)-1])]<1,Δ[z]-(Sqrt[((Δ[z])^2)-1]),Δ[z]+(Sqrt[((Δ[z])^2)-1])];
ρm[z_]=If[Abs[Δ[z]+(Sqrt[((Δ[z])^2)-1])]<1,Δ[z]+(Sqrt[((Δ[z])^2)-1]),Δ[z]-(Sqrt[((Δ[z])^2)-1])];
mp[z_]=((ρp[z]-M[z][[1,1]])/(M[z][[1,2]]));
mm[z_]=((ρm[z]-M[z][[1,1]])/(M[z][[1,2]]));

gr1=Plot[{mp[z],mm[z]},{z,-20,14.59},PlotPoints->200];
gr2=Plot[{mp[z],mm[z]},{z,16.35,46.3},PlotPoints->200,PlotRange->{-25,25}];
gr3=Plot[{mp[z],mm[z]},{z,60.6,73.75},PlotPoints->200,PlotRange->{-25,25}];
Show[gr1,gr2,gr3,PlotRange->All,AxesLabel->{"z","(SubPlus[m])(z), (SubMinus[m])(z)"}]

Export["ode5_4.eps",%];

M[z_]=FullSimplify[Φ[z+q,(1/2)].Φ[z-q,(1/2)]];

Δ[z_]=FullSimplify[(1/2)Tr[M[z]]]

sol=Flatten[Solve[{(Sqrt[z])== a + ϵ,(Sqrt[z-q])== a - ϵ},{a,ϵ}]]

ϵ ==(q/(4 a)) /. sol//Simplify

Delta[a_,ϵ_]=FullSimplify[(1/2)Tr[Φ[(((a-ϵ))^2),(1/2)].Φ[(((a+ϵ))^2),(1/2)]]]

Simplify[Delta[2n π,(q/( 4 2n π))]-1,n∈Integers]

FullSimplify[Delta[(2n+1) π,(q/(4(2 n+1)π))]+1,n∈Integers]

Planar dynamical systems

Off[General::spell1];

This command will generate the phase plot of a planar vector field. The first argument f is the
vector field. The second is a list of initial conditions. The third is the maximal time for which the solution will be computed. In addition there are some optional graphics arguments which are passed to the Show[] command.

(* Needs["VectorFieldPlots`"]; For versions before 7 replace VectorPlot by VectorFieldPlot *)

PhasePlot[f_,ic_,ti_,opts___]:=Block[{i,n=Length[ic],ff,ivp,sol,t1,t2,tmin,tmax,phaseplot,ps,pvf,vf,vfr,pev,a,ev,evs,optss},
ff=f /. {x->x[t],y->y[t]};
If[Head[ti]===List,{tmin,tmax}=ti,{tmin,tmax}={-ti,ti}];
pfv=VectorField /. {opts,VectorField->False};
pev=EigenVectors /. {opts,EigenVectors->False};
ps=PlotStyle /.{opts,PlotStyle->{}};
optss=DeleteCases[{opts},VectorField->_];
optss=DeleteCases[optss,EigenVectors->_];
optss=DeleteCases[optss,PlotStyle->_];
optss[[0]]=Sequence;
vf=Graphics[];
ev=Graphics[];
vfr=PlotRange /. {opts, PlotRange->{{-1,1},{-1,1}}};
If[pfv,vf=VectorPlot[f,{x,vfr[[1,1]],vfr[[1,2]]},{y,vfr[[2,1]],vfr[[2,2]]},Frame->False,Axes->True,Ticks->False,VectorPoints->8,VectorScale->{0.1,0.3}]];
If[pev,
a=Transpose[{D[f,x],D[f,y]}/. {x->0,y->0}];
evs=Eigensystem[a][[2]];
If[Im[evs]=={{0,0},{0,0}},
ev=ParametricPlot[evs t,{t,-100,100},PlotStyle->{Red,Red}]];
];
Off[NDSolve::"ndsz"];
Do[
ivp={x'[t]==ff[[1]],y'[t]==ff[[2]],x[0]==ic[[i,1]],y[0]==ic[[i,2]]};
sol=NDSolve[ivp,{x[t],y[t]},{t,tmin,tmax}];
{t1,t2}=sol[[1]][[1]][[2]][[0]][[1]][[1]];
phaseplot[i]=ParametricPlot[{x[t],y[t]} /. sol,{t,t1,t2},PlotStyle->ps]
,{i,1,n}];
On[NDSolve::"ndsz"];
Show[Table[phaseplot[i],{i,1,n}],vf,ev,optss]
];

This command will find all singular points of a vector field. In addition, it will
give the trace and the determinant at each singular point (from which the sing of the
real parts of the eigenvalues can be read off).

Linearize[f_,opts___]:=Block[{sol,i,jac,out,opt},
opt={DetTr,Matrix} /. {opts} /. {DetTr->True,Matrix->False};
jac=Simplify[(({{(∂_x) f[[1]], (∂_y) f[[1]]}, {(∂_x) f[[2]], (∂_y) f[[2]]}}))];
sol= Simplify[Solve[{f[[1]]==0,f[[2]]==0},{x,y}]];
Do[
out={{x,y} /.sol[[i]]};
If[opt[[1]],out=Append[out,Simplify[{Tr[jac /. sol[[i]] ],Det[jac /. sol[[i]] ]}]] ];
If[opt[[2]],out=Append[out,Simplify[MatrixForm[jac /. sol[[i]]]]] ];
Print[out]; ,{i,1,Length[sol]}];
];

Orbits and invariant sets

c=(1/10);PhasePlot[{((Cos[x])^2) (-y Cos[x]+c Sin[x]),Sin[x]+c y Cos[x]},{{0,1},{π,1},{-π,1}},50,PlotRange->{{-π,π},{-5,4.5}},Ticks->False,VectorField->True]

Linear equations

Here is the phase plot of a spiral source.

PhasePlot[{x+ y,-x +y},{{1,1},{-1,-1},{-1,1},{1,-1},{2,2},{-2,-2},{-2,2},{2,-2},{3,3},{-3,-3},{-3,3},{3,-3}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True]

Export["ode3_1a.eps",%];

and of a spiral sink.

PhasePlot[{-x- y,x -y},{{1,1},{-1,-1},{-1,1},{1,-1},{2,2},{-2,-2},{-2,2},{2,-2},{3,3},{-3,-3},{-3,3},{3,-3}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True]

Export["ode3_1b.eps",%];

No spiral part

MatrixForm[JordanDecomposition[(({{1, 1.5}, {0.5, 1}}))][[2]]]

PhasePlot[{x+0.5y,0.5x +y},{{1,1},{-1,-1},{-1,1},{1,-1},{2,1},{2,-1},{-2,1},{-2,-1},{1,2},{1,-2},{-1,2},{-1,-2},{3,-2},{2,-3},{-3,2},{-2,3}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

PhasePlot[{x+ y,0.5x +y},{{0.8164965809277261`,0.5773502691896258`},{0,1.5},{0,3.3},{0,5},{0,7},{0,10},{0,50},{-2,0},{-5,0},{-9,0},{-0.8164965809277261`,-0.5773502691896258`},{0,-1.5},{0,-3.3},{0,-5},{0,-7},{0,-10},{0,-50},{2,0},{5,0},{9,0}},10,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

Export["ode3_2a.eps",%];

PhasePlot[{-x- 0.5y,-0.5x -y},{{1,1},{-1,-1},{-1,1},{1,-1},{2,1},{2,-1},{-2,1},{-2,-1},{1,2},{1,-2},{-1,2},{-1,-2},{3,-2},{2,-3},{-3,2},{-2,3}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

PhasePlot[{-x- y,-0.5x -y},{{0.8164965809277261`,0.5773502691896258`},{0,1.5},{0,3.3},{0,5},{0,7},{0,10},{0,50},{-2,0},{-5,0},{-9,0},{-0.8164965809277261`,-0.5773502691896258`},{0,-1.5},{0,-3.3},{0,-5},{0,-7},{0,-10},{0,-50},{2,0},{5,0},{9,0}},10,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

Export["ode3_2b.eps",%];

And here is the phase plot of a saddle.

PhasePlot[{x+ y,x -y},{{1-(Sqrt[2]),1},{-1+(Sqrt[2]),-1},{1+(Sqrt[2]),1},{-1-(Sqrt[2]),-1},{1,1},{-1,-1},{-1,1},{1,-1},{2,2},{-2,-2},{-2,2},{2,-2},{3,3},{-3,-3},{-3,3},{3,-3}},3,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

Export["ode3_3.eps",%];

Finally a center

PhasePlot[{- 2y,x },{{0,1},{0,2},{0,3},{0,4},{0,5},{0,6}},π,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True]

Export["ode3_4.eps",%];

MatrixForm[JordanDecomposition[(({{1, -1}, {1, 3}}))][[2]]]

PhasePlot[{x-y,x +3y},{{-1,1},{0,1},{0,5},{0,12},{0,22},{0,50},{1,-1},{0,-1},{0,-5},{0,-12},{0,-22},{0,-50}},6,PlotRange->{{-5,5},{-5,5}},Ticks->False,VectorField->True,EigenVectors->True]

Export["ode3_5.eps",%];

Stability

PhasePlot[{-y+x(1-(x^2)-(y^2))+((x y)/(Sqrt[(x^2)+(y^2)])),x+y(1-(x^2)-(y^2))-((x^2)/(Sqrt[(x^2)+(y^2)]))},{{0,1},{-0.1,0},{-0.5,0},{0.5,0.3},{0,0.7},{1.2,0.6},{1.2,-0.2},{0.8,-1},{-1,-1},{0.2,1.2},{-1,1},{1.1,-0.4},{Cos[0.2],Sin[0.2]}},7,PlotRange->{{-1.5,1.5},{-1.5,1.5}},Ticks->False]

Newton's equations

Here is the potential of the mathematical pendulum.

U[x_]=1-Cos[x];Plot[U[x],{x,-2 π,2 π},Ticks->False]

Here is the  corresponding vector field.

Needs["VectorFieldPlots`"];

VectorFieldPlot[{y,-U'[x]},{x,-2 π,2 π},{y,-5,5},Frame->True,PlotPoints->10]

Here is the corresponding phase plot.

PhasePlot[{y,-U'[x]},{{0,0.2},{0,1},{-2π,0.2},{-2π,1},{2π,0.2},{2π,1},{0,2},{2π,-2},{2π,2},{-2π,-2},{-2π,2},{0,-2},{0,2.5},{0,-2.5},{0,3},{0,-3}},2π,PlotRange->{{-2π,2π},{-3,3}},Ticks->False]

Here is how it looks like with damping

η=1;PhasePlot[{y,-η y-U'[x]},{{-2π,1},{-3π/2,1},{-π,1},{0,1},{0,2.2},{2,1},{π,1},{2π,1},{2.3π,1},{-7,-1},{-6.2,-1},{-π,-1},{0,-1},{0,-2.2},{4,-1},{5,-1},{1.5,-1},{-2,-1}},4π,PlotRange->{{-2π,2π},{-3,3}},Ticks->False]

Here is another example.

U[x_]=((x^2)/2)- ((x^3)/3);Plot[U[x],{x,-1,2},Ticks->False]

PhasePlot[{y,-U'[x]},{{-1,0},{-0.7,0},{-0.5,0},{-0.3,0},{1.05,0.05},{1.5,0},{2,0}},8,PlotRange->{{-1,2.5},{-2,2}},Ticks->True]

Again with damping

η =1;PhasePlot[{y,-η y-U'[x]},{{-1.5,0},{-1,0},{-0.5,0},{-0.3,0},{1.05,0.05},{1.05,-0.05},{1.5,0},{2,0},{1,1},{1,2},{1,-1}},8,PlotRange->{{-1,2.5},{-2,2}},Ticks->True]

Ecological models for two species

Off[NDSolve::ndsz];Off[InterpolatingFunction::dmval];

The Volterra-Lotka equations.

PhasePlot[{(1-y)x,(x-1)y},{{0.3,0.3},{0.5,0.5},{0.7,0.7},{0.9,0.9},{1,1}},3.7,PlotRange->{{0,3.5},{0,3.5}},Ticks->False]

The Volterra-Lotka equations with limited growth.

α=1;λ=2;μ=1;xmax=4;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (x-1-μ y) y},{{0.4`,(1/8)},{1,(1/2)},{1,1},{1,3}},5,Ticks->False];
gr2=Plot[{1-λ x,((x-1)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}},AxesOrigin->{0,0}]

α=1;λ=1/2;μ=1;xmax=3;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (x-1-μ y) y},{{(1/λ),0.001`},{0.2`,1},{0.6`,1},{1,1},{2,1},{2,(1/2)}},12,Ticks->False];
gr2=Plot[{1-λ x,((x-1)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}},AxesOrigin->{0,0}]

α=1;λ=1/2;μ=3;xmax=3;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (x-1-μ y) y},{{(1/λ),0.001`},{0.2`,1},{0.6`,1},{1,1},{2,1},{2,(1/2)}},12,Ticks->False];
gr2=Plot[{1-λ x,((x-1)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}},AxesOrigin->{0,0}]

The case of competing species.

α=1;

PhasePlot[{(1-y)x,α(1-x)y},{{0.2,1},{0.6,1},{1,0.7},{1,0.3},{0.5,0.5},{2,2},{1,1.5},{1,2.2},{1.5,1},{2,1},{2.7,1},{1.5,0.6257825342012828`},{0.6257825342012828`,1.5}},12,Ticks->False,PlotRange->{{0,3},{0,3}}]

The case of competing species with limited growth.

λ=.;μ=.;α=.;
Linearize[{(1-y-λ x)x,α(1-x-μ y )y}]

α=1;λ=0.4;μ=0.55;xmax=4;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (1-x-μ y) y},{{-((1-μ)/(-1+λ μ)),-((1-λ)/(-1+λ μ))-0.01},{-((1-μ)/(-1+λ μ)),-((1-λ)/(-1+λ μ))+0.01},{0.2,0.6},{0.3,0.2},{1,2},{0.6`,2},{2,0.8},{4,1.2},{4,0.2}},21,Ticks->False];
gr2=Plot[{1-λ x,((1-x)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}}]

α=1;λ=1.5;μ=0.55;xmax=4;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (1-x-μ y) y},{{0.2,0.6},{0.3,0.2},{1,2},{0.6,2},{1.5,1.5},{4,1.2},{4,0.2}},21,Ticks->False];
gr2=Plot[{1-λ x,((1-x)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}}]

α=1;λ=1.5;μ=1.5;xmax=4;ymax=2;

gr1=PhasePlot[{(1-y-λ x) x,α (1-x-μ y) y},{{-((1-μ)/(-1+λ μ)),-((1-λ)/(-1+λ μ))-0.01`},{-((1-μ)/(-1+λ μ)),-((1-λ)/(-1+λ μ))+0.01`},{0.2`,0.6`},{0.3`,0.2`},{1,2},{0.6`,2},{2,0.8`},{4,1.2`},{4,0.2`}},21,Ticks->False];
gr2=Plot[{1-λ x,((1-x)/μ)},{x,0,xmax},PlotPoints->30,PlotStyle->Gray];
Show[gr1,gr2,PlotRange->{{0,xmax},{0,ymax}}]

The van der Pol equation

f[x_]=1(((x^3)/3)-x);
PhasePlot[{y- f[x],-x},{{0,0.1},{0,2.8},{0,-0.1},{0,-2.9}},9,Ticks->False,PlotRange->{{-3,3},{-3,3}}]

PhasePlot[{y,(1-(x^2))y-x},{{0,0.1},{0,2.8},{0,-0.1},{0,-2.9}},9,Ticks->False,PlotRange->{{-3,3},{-3.5,3.5}}]

SIR Model

β=5;γ=0.1;n0=100;
i0=1;s0=n0-1;sol={s[t],i[t],n0-i[t]-s[t]} /.NDSolve[{s'[t]==-β ((i[t])/n0)s[t],i'[t]==β ((i[t])/n0)s[t]-γ i[t],i[0]==i0,s[0]==s0},{s[t],i[t]},{t,0,20}][[1]];
Plot[sol,{t,0,14},PlotRange->All,PlotLabels->{"s","i","r"}]

If[((s0 β)/(n0 γ))>1,s0+i0-n0(γ/β)(1+Log[((s0 β)/(n0 γ))]),i0]

Clear[β,γ,i0,s0,n]

{r,r0+i0+s0-r}/.Solve[r0+i0+s0-r==s0 Exp[((-β)/γ)(((r-r0))/n)],r][[1]]//FullSimplify

Lienards equation

f[x_]=x(((x^2)/3)-1);

plf=Plot[f[x],{x,0,2},PlotStyle->Black,PlotRange->{-3,2}]

sol=NDSolve[{x'[t]==y[t]-f[x[t]],y'[t]==-x[t],x[0]==0,y[0]==2},{x[t],y[t]},{t,0,5}];
plo=ParametricPlot[{x[t],y[t]} /. sol,{t,0,3.29}];
Show[plf,plo]

Hopf bifurcation

f[x]=((x^3)/3)-μ x;

μ=-1;
PhasePlot[{y- f[x],-x},{{0,0.1},{0,2.8},{0,-0.1},{0,-2.9}},9,Ticks->False,PlotRange->{{-3,3},{-3,3}}]

μ=0;
PhasePlot[{y- f[x],-x},{{0,2.8},{0,-2.9}},100,Ticks->False,PlotRange->{{-3,3},{-3,3}}]

μ=1;
PhasePlot[{y- f[x],-x},{{0,0.1},{0,2.8},{0,-0.1},{0,-2.9}},9,Ticks->False,PlotRange->{{-3,3},{-3,3}}]

Poincare-Bendixon

Linearize[{y+(x^2)-α x(y-1+2(x^2)),-2(1+y)x}]

H[x_,y_]=(x^2)(1+y)+((y^2)/2);

{(∂_x) H[x,y],(∂_y) H[x,y]}.{y+(x^2)-α x(y-1+2(x^2)),-2(1+y)x}//Simplify

Solve[H[x,y]==1/2,y]

ContourPlot[H[x,y],{x,-1,1},{y,-1.1,1},Contours->{0.01,0.05,0.1,0.2,0.3,0.4,0.5}]

sp=Show[ContourPlot[{y==-1,y==1-2(x^2)},{x,-1,1},{y,-1.05,1},ContourStyle->Red],Graphics[{PointSize[0.015],Point[{0,0}],Point[{-1,-1}],Point[{1,-1}]}]];

α=1/4;
Show[sp,PhasePlot[{y+(x^2)-α x(y-1+2(x^2)),-2(1+y)x},{{0,-0.2}},{-20,24},PlotRange->{{-1,1},{-1,1}}]]

U[x_]=(x^2)((x^2)-1);

x[t_] =(x /.Solve[t==Simplify[(1/(Sqrt[2]))(∫_x^1)(1/(r(Sqrt[1-(r^2)])))\[DifferentialD]r,0<x<1],x][[1]]);
y[t_]=FullSimplify[D[x[t],t]];
sep=ParametricPlot[{{x[t],y[t]},{-x[t],y[t]}},{t,-10,10},PlotRange->All,PlotStyle->{Red}];
Print[x[t]];
Clear[x,y];

pp=PhasePlot[{y,-10(((((y^2)/2)+U[x]))^2)y-U'[x]},{{-1.2,0},{-0.3,0},{0.3,0}},{0,15},PlotRange->{{-1.5,1.5},{-1,1}},Ticks->False];

Show[pp,sep,Graphics[{PointSize[0.015],Point[{0,0}],Point[{(1/(Sqrt[2])),0}],Point[{((-1)/(Sqrt[2])),0}]}]]

Duffing equation

δ=0.25;
PhasePlot[{y,-δ y+x-(x^3)},{{0.01,0.01},{-0.01,-0.01}},{-17,20},Ticks->False,PlotRange->{{-2.5,2.5},{-2.4,2.4}}]

δ=0;
PhasePlot[{y,-δ y+x-(x^3)},{{2.5,0},{2,0},{1.75,0},{1.55,0},{(Sqrt[2]),0},{1.2,0},{1.35,0},{-1.35,0},{-1.2,0},{-(Sqrt[2]),0}},{-17,20},Ticks->False,PlotRange->{{-2.5,2.5},{-2.4,2.4}}]

FullSimplify[DSolve[{x'[t]==x[t](Sqrt[1-(((x[t])^2)/2)]),x[0]==(Sqrt[2])},x[t],t]]

Stable Manifold

For the PhasePlot command see section "Planar dynamical systems" above.

pp1=PhasePlot[{(1-y)x,(1-x)y+(((x-1))^2)},{{0.2,1},{0.5,1},{1,0.1},{1,1},{1,1.5},{1.6,1},{2.3,1},{1,2},{1,2.5}},12,Ticks->False,PlotRange->{{0,3},{0,3}}];

pp2=PhasePlot[{(1-y)x,(1-x)y+(((x-1))^2)},{{1.02,0.98},{1.02,1.02}},12,Ticks->False,PlotRange->{{0,3},{0,3}},PlotStyle->{Thick}];

Show[pp1,pp2,Graphics[{{Dashed,Line[{{0,0},{3,3}}]},{Dashed,Line[{{2,0},{0,2}}]}}]]

Linearize[{(1-y)x,(1-x)y},Matrix->True]

PhasePlot[{-(y-1),-(x-1)},{{0.2,1},{0.5,1},{1,0.5},{1,1},{1,1.5},{1.5,1},{2.3,1},{1,2},{1,2.5},{1.05,1.05},{0.95,1.05},{1.05,0.95},{0.95,0.95}},12,Ticks->False,PlotRange->{{0,2},{0,2}}]

pp1=PhasePlot[{(1-y)x,(1-x)y+(((x-1))^2)},{{0.2,1},{0.5,1},{1,0.1},{1,1},{1,1.5},{1.6,1},{2.3,1},{1,2},{1,2.5},{1.02,0.98},{1.02,1.02}},12,Ticks->False,PlotRange->{{0,2},{0,2}}]

pp1=PhasePlot[{-x+y+3 (y^2),y},{{0,4.4},{0,3.2},{0,2},{0,1},{2,0.3},{4,0.5},{2,-0.5},{4,-0.7},{0,-1.2},{0,-2.5},{0,-4}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False];

pp2=PhasePlot[{-x+y+3 (y^2),y},{{3/2,1},{1/2,-1},{-5,0},{5,0}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,Axes->False,PlotStyle->{Thick}];

Show[pp1,pp2,Graphics[{{Dashed,Line[{{-5,0},{5,0}}]},{Dashed,Line[{{-2.5,-5},{2.5,5}}]}}]]

PhasePlot[{-x+ y,y},{{1,1},{-1,-1},{-1,1},{1,-1},{2,2},{-2,-2},{-2,2},{2,-2},{3,3},{-3,-3},{-3,3},{3,-3}},5,PlotRange->{{-5,5},{-5,5}},Ticks->False,EigenVectors->True,Axes->False]

FPU experiment

n=32;
m=1;k=1;α=1/6;τ=2(n+1);
tmax=500τ;
U[x_]=(k/2)((x^2)+α (x^3));
ode=Simplify[Flatten[Table[{(∂_t) p[j,t] ==U'[q[j+1,t]-q[j,t]]- U'[q[j,t]-q[j-1,t]],(∂_t) q[j,t] == ((p[j,t])/m)},{j,n}] /.{q[0,t]->0,q[n+1,t]->0}]];
vars=Flatten[Transpose[Table[{p[j,t],q[j,t]},{j,1,n}]]];

Linear modes; Initial condition given by the first linear mode.

v[j_]=Table[(Sqrt[(2/(n+1))])Sin[((j π)/(n+1))i],{i,n}];
init=Thread[(vars/. t->0)==  Flatten[{Table[0,{i,n}],v[1]}]];

Solve the ODE

eqs=Flatten[{ode,init}];sol=NDSolve[eqs,vars,{t,0,tmax}][[1]];

Here is a plot for the first particle over one period of the linear system.

Plot[{p[1,t],q[1,t]}/.sol,{t,0,τ}]

Decompose the solution into the linear modes:

a[j_,t_]=Table[q[i,t],{i,n}].v[j]/.sol; b[j_,t_]=Table[p[i,t],{i,n}].v[j]/.sol;

Evaluate the (linear) energy corresponding to the first four modes at multiples of the linear period:

ListPlot[Table[Table[{(t/τ),(1/(2 m))((b[j,t])^2)+k(1- Cos[(( π j)/(n+1))])((a[j,t] )^2)}/.sol,{t,0,tmax,τ}],{j,4}]]

The Lorenz system

tmax=20;σ=10;ρ=28;β=(8/3);orbit=NDSolve[{x'[t]==-σ (x[t]-y[t]),y'[t]==-x[t] z[t]+ρ x[t]-y[t],z'[t]==x[t] y[t]-β z[t],x[0]==30,y[0]==10,z[0]==40},{x,y,z},{t,0,tmax},MaxSteps->5000];
gr=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.orbit],{t,0,tmax},PlotPoints->2000,Axes->False,PlotRange->All]

Show[gr,Graphics3D[{RGBColor[0,0,1],PointSize[0.01`],Point[{-(Sqrt[β (ρ-1)]),-(Sqrt[β (ρ-1)]),ρ-1}],Point[{(Sqrt[β (ρ-1)]),(Sqrt[β (ρ-1)]),ρ-1}]}]]

Clear[σ,r,b]

f[x_,y_,z_]:={-σ (x-y),r x-x z-y,x y-b z};

Solve[Thread[f[x,y,z]=={0,0,0}],{x,y,z}]

jac=Simplify[(({{(∂_x) f[x,y,z][[1]], (∂_y) f[x,y,z][[1]], (∂_z) f[x,y,z][[1]]}, {(∂_x) f[x,y,z][[2]], (∂_y) f[x,y,z][[2]], (∂_z) f[x,y,z][[2]]}, {(∂_x) f[x,y,z][[3]], (∂_y) f[x,y,z][[3]], (∂_z) f[x,y,z][[3]]}}))] ;

MatrixForm[jac]

MatrixForm[jac/. {x->0,y->0,z->0}]

Eigenvalues[jac /. {x->0,y->0,z->0}]//FullSimplify

MatrixForm[jac/. {x->(Sqrt[b(r-1)]),y->(Sqrt[b(r-1)]),z->r-1}]

Apply the Routh-Hurwitz test

a=Det[z IdentityMatrix[3]-%]//Drop[Reverse[CoefficientList[#,z]],1]&

{a[[1]]>0,Det[(({{a[[1]], 1}, {a[[3]], a[[2]]}}))]>0,Det[(({{a[[1]], 1, 0}, {a[[3]], a[[2]], 1}, {0, 0, a[[3]]}}))]>0}//Simplify

FullSimplify[1+b+σ>0&&b (r (1+b-σ)+σ (3+b+σ))>0&&(b^2) (-1+r) σ (r (1+b-σ)+σ (3+b+σ))>0,Assumptions->{σ>0,b>0,r>1}]

L[x_,y_,z_]=α (x^2)+β (y^2)+γ (z^2);

D[L[x[t],y[t],z[t]],t] /. {x'[t]->-σ (x[t]-y[t]),y'[t]->r x[t]-x [t]z[t]-y[t],z'[t]->x[t] y[t]-b z[t]}//FullSimplify

-2 σ (r (((x[t]-y[t]))^2)+(1-r)((y[t])^2)+b ((z[t])^2))==% /. {γ->β,α->r,β->σ} //Simplify

σ = 0

Clear[ρ,β]

{λ1,λ2}=Eigenvalues[(({{-1, -x0}, {x0, -β}}))]

U=(({{λ1+β, λ2+β}, {x0, x0}})); Inverse[U].(({{-1, -x0}, {x0, -β}})).U//MatrixForm//FullSimplify

Clear[λ1,λ2]

(P[t_]=FullSimplify[U.(({{Exp[λ1 t], 0}, {0, Exp[λ2 t]}})).Inverse[U]])//MatrixForm

(∫_0^t)P[t-s].{ r x0, 0}\[DifferentialD]s//Simplify

The logistic map

Clear[μ]

ShowWeb[f_,xstart_,nmax_,opts___]:=Module[{x,xmin,xmax,delta,graph,web,lines},
lines=Type /. {opts,Type->Line};
x[0]:=xstart;
x[n_]:=x[n]=f[x[n-1]];
If[x[0]==x[1],Return[Print["Fixed Point!"]]];
web=Flatten[Table[{{x[n],x[n]},{x[n],x[n+1]}},{n,0,nmax}],1];
xmax=Max[web];xmin=Min[web];delta=0.1(xmax-xmin);
graph=Plot[{f[x],x},{x,xmin-delta,xmax+delta}];
Show[graph,Graphics[lines[web]]]
] /;(IntegerQ[nmax]&&Positive[nmax])

Arrows[l_List]:={Arrowheads[Small],Table[Arrow[{l[[i]],l[[i+1]]}],{i,Length[l]-1}]};

ShowWeb[1#(1-#)&,0.4,20]

ShowWeb[2#(1-#)&,0.2,20]

ShowWeb[2#(1-#)&,0.2,4,Type->Arrows]

ShowWeb[3.1#(1-#)&,0.4,20]

Solve[μ x(1-x)==x,x]

Solve[-(-1+x) x (μ^2) (1-x μ+(x^2) μ)==x,x]//FullSimplify

BifurcationList[f_,x0_,{μ_,μ0_,μ1_},opts___]:=
Block[{Nmin,Nmax,Steps,Delta,mychop},
{Nmin,Nmax,Steps,Delta}={Nmin,Nmax,Steps,Delta} /. {opts} /. {Nmin->300,Nmax->500, Steps->300,Delta->0.002};
If[Delta>0,mychop=Compile[{x},  Evaluate[Round[x /Delta]Delta]],mychop[x_]=Chop[x]];
Flatten[Table[
Module[{x,l},
x=Nest[f,x0,Nmin];
l=NestList[f,x,Nmax-Nmin];
l=Map[mychop,l];
l=Union[l];
Map[{μ,#}&,l] ],
{μ,μ0,μ1,((μ1-μ0)/Steps)}],1]
];

list=BifurcationList[μ#(1-#)&,0.4,{μ,2.95,4.}];

Length[list]

StripBifurcationList[list_,del_]:=Block[{d,l1,l2},
l1=Transpose[list][[1]];
l2=Transpose[list][[2]];
d=((Max[l2]-Min[l2])/del);
Union[Transpose[{l1,Map[d Round[#/d]&,l2]}]]
]

Length[list=StripBifurcationList[list,300]]

ListPlot[list,PlotStyle->{PointSize[0.002`]},PlotRange->All,Axes->False,AspectRatio->0.8`]

μ=.

f[x_]:=μ x (1-x);p=1-(1/μ);

g[x_]=Nest[f,x,3]

g[p]

eqn=Simplify[g[x] ==x/. μ->4]

sol=Simplify[Union[Solve[eqn,x]]]

FullSimplify[ComplexExpand[(1/24) (14+((2 (7^(2/3)))/(((-(1/2)+((3 I (Sqrt[3]))/2)))^(1/3)))+(2^(2/3)) (((-7+21 I (Sqrt[3])))^(1/3)))]]

FullSimplify[((((1/2) (1+I (Sqrt[3]))))^(1/3))]

FullSimplify[(1/8) (4+(2/(((-1))^(1/9)))+2 (((-1))^(1/9)))]

FullSimplify[ComplexExpand[%]]

Simplify[f[%] /. μ->4]

Simplify[f[%] /. μ->4]

N[%]

Chop[N[%]]

f[f[0.6112604669781573`]]

eqs={x==g[x],g'[x]==1}

sol=Union[Solve[eqs,{x,μ}]];

sol[[8]]

Chop[N[%]]

FullSimplify[((((1/2) (25-22 (Sqrt[2])+I (Sqrt[43011-29700 (Sqrt[2])]))))^(1/3))]

FullSimplify[(1/(2^(1/3)))(((-1+I 3 (Sqrt[3])))^(1/3))(-1+2 (Sqrt[2]))]

Simplify[(1/21) (10+(Sqrt[2]))-(((-1+2 (Sqrt[2])) (((-(1/2)+((3 I (Sqrt[3]))/2)))^(1/3)))/(3 (7^(2/3))))+((-9+4 (Sqrt[2]))/(3 ((7 )^(1/3))(-1+2 (Sqrt[2])) (((-(1/2)+((3 I (Sqrt[3]))/2)))^(1/3))))]

Simplify[ComplexExpand[%]]

FullSimplify[%]

N[%]

sol=Union[Solve[{y==f[x],z==f[y],x==f[z],f'[x]f'[y]f'[z]==1},{x,y,z,μ}]];

Chop[N[sol[[6]]]]

FullSimplify[sol[[6]]]

Solve[f[x]==x,x]

Simplify[((f[x]-p)/(x-p))]

FullSimplify[((f[f[x]]-p)/(x-p))]

sol=Solve[f[f[x]]==x,x]//FullSimplify

D[f[f[x]],x] /.sol//Simplify

Solve[4+2 μ-(μ^2)==-1]

f[x_]=μ x (1-x);

g[y_]=x /.Solve[f[x]==y,x][[1]];

g[1] /. μ->4

The tent map

μ=0.9;
F[x_]=μ(1-2Abs[x-(1/2)]);
F[n_,x_]:=Nest[F,x,n]

Plot[Evaluate[Table[F[n,x],{n,3}]],{x,0,1},PlotRange->{0,1},PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]

Homoclinic orbit

μ=5;
x0=Nest[((1/2)- (Sqrt[(1/4)- ((#)/μ)]))&,1.,5];
ShowWeb[μ#(1-#)&,x0,6]

Melnikov method

Duffing equation

q[t_]=(Sqrt[2])Sech[t];p[t_]=D[q[t],t];

Simplify[(∫_(-∞)^∞)p[s](δ p[s] + γ Cos[ω(s-t)])\[DifferentialD]s, ω∈Reals]

Mathematical Pendulum

q[t_]=4ArcTan[Tanh[(t/2)]] ;
p[t_]=FullSimplify[D[q[t],t]];

(∫_(-∞)^∞)p[s]Cos[s-t]\[DifferentialD]s