This notebook shows the relationship between the Sturm-Liouville Eigenvalue problem
(d/dx)(p(x) (du/dx)) + q(x) u + λ σ(x) u = 0  plus boundary conditions
for eigenfunctions u(x) and the generalized eigenvalue problem Av = λBv for eigenvectors v. 
While the generalized eigenvalue problem is useful for understanding why the eigenvalues are real, and the eigenvectors (of eigenfunctions) are orthogonal,  the Mathematica calculation is much faster if we find the eigensystem of (B^(-1))Av = λv rather than Av = λBv.

The companion notebook (SturmLiouville.nb) does the faster calculation.
2019-12-06 version

Clear[p0, p1, p2, p3, p4,p5];
σMat = DiagonalMatrix[{σ1, σ2, σ3, σ4}];
Print["σMat = ", MatrixForm[σMat], " evaluated at (Subscript[(OverBar[x]), 1]), (Subscript[(OverBar[x]), 2]), (Subscript[(OverBar[x]), 3]), (Subscript[(OverBar[x]), 4]), midpoints of n = 4 subintervals.  qMat is similar."]
pMat = DiagonalMatrix[{p0,p1, p2,p3, p4}];
Print["pMat = ", MatrixForm[pMat], "evaluated at (Subscript[x, 0]) = a, (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4]) = b, endpoints of n = 4 subintervals."]
Dtall = {{α, 0, 0, 0},{-1, 1, 0, 0},{0,-1,1,0},{0,0,-1,1}, {0,0,0,β}};
Dwide = {{-1,1, 0, 0, 0},{0, -1,1,0, 0}, {0,0,-1,1,  0}, {0,0,0,-1,1}};
Print["Dwide = ", MatrixForm[Dwide], " approximates derivatives at midpoints given function at endpoints."]
Print["Dtall = ",MatrixForm[Dtall], "  α and β approximates derivatives at x = a and x = b, depending on boundary conditions."]
Print["(d/dx) p(x)(d/dx) = (1/((Δx)^2)) Dwide.pMat.Dtall = (1/((Δx)^2)) ",MatrixForm[Dwide.pMat.Dtall]]

sturmLiouville := Module[{},
Δ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 *)
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[{Dwide.pMat.Dtall/((Δx)^2) + qMat, -σMat}, -numE];
Print["Generalized evec calculation took ", time, " seconds with n = ", n];
evals = Reverse[evals]; 
Print["The eigenvalues are ", evals];
p1 = ListPlot[evecs[[1]]*Sign[evecs[[1,1]] ] ] ;
{time,{evals, evecs}} = Timing@Eigensystem[-σInvMat.(Dwide.pMat.Dtall/((Δx)^2) + qMat), -numE];
Print["Evec calculation using σInv took ", time, " seconds with n = ", n];
 (* We want to list them from smallest to largest *)
evals = Reverse[evals]; 
evecs = Reverse[evecs];
Print["The eigenvalues are ", evals];
p1 = plotEvec[evecs[[1]], PlotStyle -> Black  ] ;
p2 = plotEvec[evecs[[2]]  ] ;
p3= plotEvec[evecs[[3]], PlotStyle -> Red  ] ;
p4 = plotEvec[evecs[[4]], PlotStyle -> Blue  ] ;
out1 = Show[p1, p2, p3, p4,PlotRange -> All];
out2 = GraphicsGrid[{{p1, p2}, {p3, p4}}, ImageSize->800];
out1
]
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]
];

n = 500; (* number of subivisions of [a, b] *)
(* With Bessels equation with m = 0, 1, or 2, with n = 1000 was OK but n = 2000 took a long time 
Actually it is slow even for the simplest.  Does the generalized eigenvector slow it down? *)
numE = 4; (* number of eiegvalues/eigenvectors computed *)
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;
(* u'' + lambda u = 0.  The eigenvalues are n^2 for positive integers n. *)
p[x_] = 1;  
q[x_] =0; 
σ[x_] = 1;
(*Data is above here *)
sturmLiouville

(* pg. 338 (7.10.11) for Legendre Polynomials d/dx (Sin[x]) dg/dx) -m^2/Sin[x] g + λ Sin[x] g = 0  *)
p[x_] = Sin[x];  
q[x_] =-0^2/Sin[x];  (* q(x) = -m^2/Sin[x]  *)
σ[x_] = Sin[x];
a = 0; (* Interval is [a, b], or a ≤ x ≤ b *)
b = π;
(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 = 500; (* number of subivisions of [a, b] *)
(* With m = 1000 and m = 0, the Legendre polynomial calculation took 15.1 seconds with the regular eigensystem calculation and 0.016 with the generalized eigensystem calculation. The saved result is below*)
numE =6; (* number of eiegvalues/eigenvectors computed.  It gets rounded up to a multiple of 3. *)
sturmLiouville

This result is saved, showing result for n = 1000.  I don’t want to wait 15 seconds for this.