% Solves mixed problem u''(x) = f(x), 0 < x < 1
% u(0) = 0, u'(1) = sigma
f = @(x) pi^2*sin(pi*x);
sigma = -pi;
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
b(n) = b(n) + sigma;
A(1,:) = []; % remove boundary conditions
A(:,1) = [];
b(1) = [];
u = A \ b; % solve system
u = [0;u];
% plot
plot(x,u)