Sturm-Liouville
This notebook estimates eigenvalues and eigenfunctions of a Sturm-Liouville problem by approximating the operator with boundary conditions as a symmetric matrix.
version 2020-02-04

sturmLiouville := Module[{imageSize = 600},
Δx = (b-a)/n;
xList =Table[a + (b-a)(i-.5)/n,{i, n}]; (* cell centers *)
pList = Table[p[a + (b-a)(i-1)/n], {i, n+1}] ;(* evaluated at endpoints of subintervals *)
DtTL = ((2 (Subscript[β, 1]) Δx)/((Subscript[β, 1]) Δx - 2(Subscript[β, 2]) ));  (* Dtall Top Left element *)
DtBR =((-2 (Subscript[β, 3]) Δx)/((Subscript[β, 3]) Δx + 2(Subscript[β, 4]) ));  (* Dtall Bottom Right element *)
numE = Max[numE, plotsPerRow*rowsOfPlots];
qList = Table[q[a + (b-a)(i-.5)/n],{i, n}]; (* q and σ are evaluated at midpoints of the interval *)
σList = Table[σ[a + (b-a)(i-.5)/n],{i, n}];
σInvList = Table[1/σ[a + (b-a)(i-.5)/n],{i, n}];
Dwide = N@SparseArray[{Band[{1,1}]->-1, Band[{1,2}]-> 1},{n,n+1}];
Dtall = N@SparseArray[{{1,1}-> DtTL,{n+1,n}-> DtBR,Band[{1,1}]->1, Band[{2,1}]-> -1},{n+1,n}];
pMat = N@SparseArray[Band[{1,1}]-> pList, {n+1, n+1}];
qMat = N@SparseArray[Band[{1,1}]-> qList, {n, n}];
σMat = N@SparseArray[Band[{1,1}]-> σList, {n, n}]; 
σInvMat = N@SparseArray[Band[{1,1}]-> σInvList, {n, n}]; (* The N[ ] makes it a floating point matrix *)
{time,{evals, evecs}} = Timing@Eigensystem[-σInvMat.(Dwide.pMat.Dtall/((Δx)^2) + qMat), -numE];
Print["Eigensystem calculation took ", time, " seconds with n = ", n];
 (* We want to list them from smallest to largest, in absolute value *)
evals = Reverse[evals]; 
evecs = Reverse[evecs];
Print["The eigenvalues are ", evals];
pTable= Table[plotEvec[evecs[[i]] ], {i, 1, rowsOfPlots*plotsPerRow}];
GraphicsGrid[Partition[pTable, plotsPerRow], ImageSize -> imageSize]
]
plotEvec[list_, options___ ] := Module[{plotPts = 400},
step = Max[1, Floor[n/plotPts] ];
sign = 1;
If[list[[step]] < 0, sign= -1];(* this makes first component plotted positive  *)
xyList = Transpose[{xList, sign*list}];
xyListToPlot = xyList[[ ;; ;; step ]];
AppendTo[xyListToPlot, xyList[[-1]] ];
ListPlot[xyListToPlot, Joined -> True, options]
];

(* ϕ'' + λ ϕ = 0, ϕ(0) = ϕ(ℼ) = 0 *)
p[x_] = 1;  
q[x_] =0;  
σ[x_] = 1;
a = 0; (* Interval is [a, b], or a ≤ x ≤ b *)
b = Pi;
(Subscript[β, 1])=1;  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=0;
(Subscript[β, 3])=1; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=0;  (* (Subscript[β, 3]) (Subscript[β, 4]) > 0 gives Newton's law of cooling, (Subscript[β, 3]) (Subscript[β, 4]) < 0 gives "unphysical case" *)
n = 1000; (* number of subivisions of [a, b] *)
numE =6; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 2; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 2; (* plots this may rows *)
sturmLiouville

The eigenvalues of ϕ’’ + λ ϕ = 0, ϕ(0) = 0, ϕ(ℼ) + ϕ’(ℼ) = 0 are approximately 0.6, 2.8, 6.8, etc.  They are the zeros of the function graphed here:

Plot[Sin[(Sqrt[λ]) Pi] + (Sqrt[λ]) Cos[(Sqrt[λ]) Pi], {λ, 0, 7}]

FindRoot[Sin[(Sqrt[λ]) Pi] + (Sqrt[λ]) Cos[(Sqrt[λ]) Pi] == 0, {λ, 7}]

The following cell computes these eigenfunctions numerically.

(* ϕ'' + λ ϕ = 0, , ϕ(0) = 0,  ϕ(ℼ)+ ϕ'(ℼ) = 0 *)
p[x_] = 1;  
q[x_] =0;  
σ[x_] = 1;
a = 0; (* Interval is [a, b], or a ≤ x ≤ b *)
b = Pi;
(Subscript[β, 1])=1;  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=0;
(Subscript[β, 3])=1; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=1;  (* (Subscript[β, 3]) (Subscript[β, 4]) > 0 gives Newton's law of cooling, (Subscript[β, 3]) (Subscript[β, 4]) < 0 gives "unphysical case" *)
n = 1000; (* number of subivisions of [a, b] *)
numE =6; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 2; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 2; (* plots this may rows *)
sturmLiouville

(* ϕ'' + λ ϕ = 0 *)
p[x_] = 1;  
q[x_] =0;  
σ[x_] = 1;
a = 0; (* Interval is [a, b], or a ≤ x ≤ b *)
b = 1;
(* only ratio of βs matters. We can have just 1 parameter, with 0 ≤ θ < π. *)
(Subscript[θ, a ])= 3π/8; (* θ = 0 is 0-Dirichlet, θ = π/2 is 0-Neumann, θ > π/2 gives "unphysical cases" *)
(Subscript[β, 1])=Cos[(Subscript[θ, a])];  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=-Sin[(Subscript[θ, a])];
(Subscript[θ, b ])=3π/8;  (* (Subscript[θ, a ]) = (Subscript[θ, b ])gives symmetry, and eigenfunctions are either even or odd *)
(Subscript[β, 3])=Cos[(Subscript[θ, b])]; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=Sin[(Subscript[θ, b])];  
n = 500; (* number of subivisions of [a, b] *)
numE =6; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 2; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 2; (* plots this may rows *)
sturmLiouville

(* pg. 302 Haberman PDE text 5E,(section 7.7.7) d/dr (r dϕ/dr) -m^2/r ϕ + λ r ϕ = 0  
Bessel Functions have p(x) = x, q(x) = -m^2/x, σ(x) = x*)
p[x_] = x;  (* Bessel Functions *)
q[x_] =-9/x;  (* q(x) = -m^2/x for sin(m θ)  *)
σ[x_] = x;
a = 0; (* Interval is [a, b], or a ≤ x ≤ b *)
b = 1;
(Subscript[β, 1])=1;  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=0;
(Subscript[β, 3])=1; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=0;
n = 1000; (* number of subivisions of [a, b] *)
numE =6; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 3; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 2; (* plots this may rows *)
sturmLiouville

(* pg. 332 in Haberman 5E,  (section 7.10.3) for Legendre Polynomials d/dx ((1 - x^2) dϕ/dx) -m^2/(1-x^2) ϕ + λ ϕ = 0  *)
p[x_] = 1-x^2;  
m = 1;
q[x_] =-m^2/(1-x^2);  (* q(x) = -m^2/(1-x^2), m = 0 gives Legendre Polynomials, other integers give associated Legendere functions, which are polynomials for m even. *)
σ[x_] = 1;
a = -1; (* Interval is [a, b], or a ≤ x ≤ b *)
b = 1;
(Subscript[β, 1])=1;  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=0;
(Subscript[β, 3])=1; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=0;
n = 1000; (* number of subivisions of [a, b] *)
numE =6; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 3; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 2; (* plots this may rows *)
sturmLiouville

(* ϕ'' - λ x ϕ = 0, eigenvalue problem that comes from the Tricomi equation*)
p[x_] = 1;  
q[x_] =0;  
σ[x_] = -x;
a = -2; (* Interval is [a, b], or a ≤ x ≤ b *)
b = 1;
(Subscript[β, 1])=1;  (* (Subscript[β, 1]) ϕ(a) + (Subscript[β, 2]) ϕ'(a) = 0 *)
(Subscript[β, 2])=0;
(Subscript[β, 3])=1; (* (Subscript[β, 3]) ϕ(b) + (Subscript[β, 4]) ϕ'(b) = 0 *)
(Subscript[β, 4])=0;  
n = 1000; (* number of subivisions of [a, b] *)
numE =9; (* number of eigevalues/eigenvectors computed. *)
plotsPerRow = 3; (* Plots this many eigenfunctions in each row*)
rowsOfPlots = 3; (* plots this may rows *)
sturmLiouville
Chop@MatrixForm[evecs.σMat.Transpose[evecs]] (* eigenfunctions are orthogonal with the weight matrix σMat *)