function [] = linear_fe(N)

% randomly select N nodes in the interval from 0 to1
x = rand(N,1);

% sort these nodes 
x = sort(x);

% scale so that x(1) = 0, x(end) = 1;
x = x - min(x);
x = x/max(x);

h = x(2:end) - x(1:end-1);

% Initialize stiffness matrix using for loops
K = zeros(N-1,N-1);

% first row:
K(1,1) = 1/h(1) + 1/h(2);
K(1,2) = - 1/h(2);

% last row
K(end,end) = 1/h(end);
K(end,end-1) = -1/h(end);

% remaining rows
for i = 2:(N-2)
    K(i,i-1) = - 1/h(i);
    K(i,i) = 1/h(i) + 1/h(i+1);
    K(i,i+1) = -1/h(i+1);
end


% Can initialize stiffness matrix using matlab vector notation
d1 = [-1./h(2:end); 0];
d2 = [ (1./h(1:end-1) + 1./h(2:end)); 1./h(end)];
d3 = -1./h;
K2 = spdiags([d1,d2,d3],-1:1, N-1, N-1);

% difference between the two:
max(max(abs(K - K2)))

% right hand side:
f = sin(2*pi*x(2:end));

F = [ 1/2*(h(1:end-1)+h(2:end)).*f(1:end-1);
      h(end)/2*f(end)];

% coefficients:
U = K\F;

% plot solution.
figure(1); clf
plot(x,sin(2*pi*x)/(2*pi)^2-x/2/pi,'r-',...
     x(2:end), U, 'bo-');
legend('exact','linear fe');
