% This function solves the Sturm-Liouville problem using FEM.
% Trapetzoidal numerical integration rule is used to compute the
% stiffness matrix and load vector.
% This function assumes homogenous Dirichlet boundary conditions.
function [x, v] = sl_fem(afn, ffn, h)
npts = length(h) + 1;
x = 0;
for k=1:length(h)
  x(k+1) = x(k) + h(k);
end

av = eval([afn '(x)']);

d = (av(1:npts-2) + av(2:npts-1)) ./ (2*h(1:npts-2)) + ... 
	 (av(3:npts) + av(2:npts-1)) ./ (2*h(2:npts-1));
ud = -(av(3:npts-1) + av(2:npts-2)) ./ (2*h(2:npts-2));

A = spdiags([[ud, 0]', d', [0, ud]'], [-1 0 1], npts-2, npts-2);
b = (h(1:npts-2) + h(2:npts-1)) .* eval([ffn '(x(2:npts-1))']) / 2;

v = [0; A \ b'; 0];


