f = @(x) pi^2*sin(pi*x);
n = 100; % number of nodes
% create grid
x = linspace(0,1,n); 
% vector of element sizes
h = diff(x);
% initialize load and stiffness
A = zeros(n); 
b = zeros(n,1); 
% element loop
for el=1:n-1 
k = el:el+1;
A(k,k) = A(k,k) + [1,-1;-1,1]/h(el);
% use 1-point 'barycentric' (midpoint) quadrature
b(k) = b(k) + f(mean(x(k)))*h(el)/2;
end
A([1,n],:) = []; % remove boundary conditions
A(:,[1,n]) = [];
b([1,n]) = [];
u = A \ b; % solve system
u = [0;u;0];
% plot
plot(x,u)