% regarding the Sturm-Liouville type equation: 
% (-p(x)u')' + q(x)u = f(x), 0 < x < 1, u(0)=0
% values below correspond to Legendre equation
p = @(x) x.^2-1;
q = @(x) 1;
f = @(x) 0;
% number of nodes
n = 100; 
% create grid
x = linspace(0,1,n); 
% vector of element sizes
h = diff(x);
% initialize load, mass and stiffness
A = zeros(n); 
M = zeros(n);
b = zeros(n,1); 
% element loop - this time need quadrature for stiffness and mass.
for el=1:n-1 
k = el:el+1;
% use 1-point 'barycentric' (midpoint) quadrature
bary = mean(x(k));
A(k,k) = A(k,k) + p(bary)*[1,-1;-1,1]/h(el);
% exercise: check that this is the correct expression for the mass matrix!
M(k,k) = M(k,k) + q(bary)*[2,1;1,2]*h(el)/6;
b(k) = b(k) + f(bary)*h(el)/2;
end
% the above code could be used to solve, e.g. Dirichlet problem by
% inputting two-sided boundary conditions. Instead, for some fun:
% implement boundary conditions on one-side only:
A(1,:) = []; 
A(:,1) = [];
M(1,:) = [];
M(:,1) = [];
b(1) = [];
% boundary conditions only on one-side, because we now solve the eigenvalue
% problem Au = lambda.Mu, here implemented as Mu = D.Au to pick the
% smallest eigenvalues (rather than largest).
[V,D] = eigs(-M,A,5);
% the true eigenvalues are n(n+1), where n=1,3,5,..., compare:
lambda=1./diag(D)
V = [zeros(1,5);V];
% the solutions to this problem are the Legendre polynomials of odd degree:
plot(x,V(:,1)/V(end,1))
hold on
plot(x,V(:,2)/V(end,2))
plot(x,V(:,3)/V(end,3))
%plot(x,V(:,4)/V(end,4))
%plot(x,V(:,5)/V(end,5))