% P1 finite elements in 1D
% Fernando Guevara Vasquez 2010
n=5; % number of elements

x = linspace(0,1,n+1)';
xf = linspace(0,1,1000)';
f = @(x) exp(x);
utrue = @(x) -exp(x) + (exp(1)-1)*x + 1;

% local to global map
i = [1:n;2:n+1]';

% local stiffness matrix in the parent element
Kloc = [1 -1; -1 1];

% build the system matrix
K = sparse(n+1,n+1);
for e=1:n,
 K(i(e,:),i(e,:)) = K(i(e,:),i(e,:)) + Kloc/(x(e+1)-x(e));
end;

K(1,:) = 0; K(1,1)=1;
K(n+1,:) = 0; K(n+1,n+1)=1;

% local mass matrix
Floc = [1/3 1/6; 1/6 1/3];
% build right hand side
b=zeros(n+1,1);
for e=1:n,
 b(i(e,:)) = b(i(e,:)) + (x(e+1)-x(e))*Floc*f(x(i(e,:)));
end;
b(1)=0; b(n+1)=0;

% solve
u = K\b;

% plot and compare
plot(x,u,xf,utrue(xf));
xlabel('x'); ylabel('u');
