% this solves -u''(x) + u(x) = x, 0 < x < 1, u(0) = u(1) = 0
% using the finite element method with N elements of size h = 1/N
function [approx,err]=fem1(N)
% L: domain size
L=1;
% h: element size
h=L/N;
% make stiffness matrix;
% diagonal elements: (phi_i',phi_i') = 2/h, (phi_i,phi_i) = 2h/3
K=(2/h+2*h/3)*eye(N-1);
% superdiagonal: (phi_i',phi_i+1') = -1/h, (phi_i,phi_i+1) = h/6
K=K+(h/6-1/h)*diag(ones(N-2,1),1);
% subdiagonal: (phi_i+1',phi_i') = -1/h, (phi_i+1,phi_i) = h/6
K=K+(h/6-1/h)*diag(ones(N-2,1),-1);
% load vector
f=h^2*(1:N-1)';
% solve the system Ku = f for the values of u at the interior nodes
% pad with zeros at either end to form the complete solution
approx=[ 0; K\f; 0 ];
figure(1)
clf
% plot solution at nodes and at several other points, to show error
% inside each element
xf=(0:h/10:1)';
exactf=xf-sinh(xf)/sinh(1);
plot(xf,exactf)
hold on
% plot approximate solution
x=h*(0:N)';
plot(x,approx,'r--')
xlabel('x')
ylabel('u(x)')
legend('exact','approx')
% estimate error using infinity vector norm on nodal values
exact=x-sinh(x)/sinh(1);
err=max(abs(exact-approx))./max(abs(exact));
