% solves inhomogeneous Dirichlet problem u''(x) = f(x), 0 < x < 1
% u(0) = a1, u(1) = a2; 
f = @(x) pi^2*cos(pi*x);
a1 = 1;
a2 = -1;
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
% construct lift
R = [a1;zeros(n-2,1);a2];
b = b - A*R;
% implement boundary conditions
A([1,n],:) = []; 
A(:,[1,n]) = [];
b([1,n]) = [];
u = A \ b; % solve homogeneous system
u = [a1;u;a2];
% plot
plot(x,u)