function [x, v] = fem(M,f)

% piecewise finite element method for
% -u_xx = f
% u(0) = u(1) = 0

% Note: h = 1/M





h = 1/M;
x = linspace(0,1,M+1);




% form the linear system

A = (2/h)*diag(ones(M-1,1));
A = A + (-1/h)*diag(ones(M-2,1),1) + (-1/h)*diag(ones(M-2,1),-1);

b = zeros(M-1,1);
for i = 1:M-1
    
    phi = @(y) (y-x(i))./h;
    b(i) = b(i) + quad(@(y) f(y).*phi(y), x(i),x(i+1));
    
    phi = @(y) (x(i+2)-y)./h;
    b(i) = b(i) + quad(@(y) f(y).*phi(y), x(i+1),x(i+2));
    
end


% solve the linear system

vtmp = A\b;


% output

v = [0; vtmp; 0];