% Create uniform partition with N nodes for (0,1)
N = 1000;
x = linspace(0,1,N);

% define the load function.
f = @(x)(ones(size(x)));

% Initialise the matrix Ahat and vector bhat. 
Ahat = sparse(N,N);
bhat = zeros(N,1);

% loop over the intervals 
for k = 1:(length(x)-1)
    
    % extract endpoints of the interval
    x1 = x(k);
    x2 = x(k+1);
    
    % evaluate length of interval k.
    len = x2-x1;

    % evaluate derivatives of basisfunctions on interval k.
    dphi(1) = 1/(x1-x2);
    dphi(2) = 1/(x2-x1);
    
    % midpoint quadrature points 
    t = (x1+x2)/2; w = x2-x1;
    
    % evaluate values of basisfunctions 
    % source term at integration points
    phi(:,1) = (t-x2)./(x1-x2);
    phi(:,2) = (t-x1)./(x2-x1);
    
    fval = f(t);
   
    % enumerate the basisfunctions on interval k.
    enum([1 2]) = [k k+1];
    for i=1:2
        % evaluate intergrals related to b.
        bhat(enum(i) ) = bhat(enum(i) ) + dot(fval.*phi(:,i),w);
        
        for j=1:2
            % evaluate integral related to A
            Ahat(enum(i),enum(j)) = Ahat(enum(i),enum(j)) + dphi(i)*dphi(j)*len;        
        end
    end    
end

% remove basisfunction 1 and N+1 from the system
A = Ahat(2:(N-1), 2:(N-1) );
b = bhat(2:(N-1),1);

u(1,1) = 0;
u(N,1) = 0;
u(2:(N-1),1) = A\b;

% plot the solution
plot(x,u);