function [tp,yp] = fem1(N,Np);
%
%   function [tp,yp] = fem1(N,Np);
%
% This function implements the finite element method with piecewise
% linear continuous basis functions for the PDE y'' + y = -t subject
% to homogeneous Dirichlet boundary conditions. N denotes the number 
% of basis functions, and Np denotes the size of the arrays for tp an
% the function values of the approximation yp for plotting.
%

t  = linspace(0,1,N+2);
tp = linspace(0,1,Np);
h  = diff(t);

% Right-hand side evaluated at the mesh points

g = -t;

% Create the stiffness matrix

A = zeros(N,N);

A(1,1) = -1.0 / h(1)   + h(1)   / 3.0;
A(N,N) = -1.0 / h(N+1) + h(N+1) / 3.0;

for i=1:N-1
    A(i  ,i  ) = A(i  ,i  ) - 1.0/h(i+1) + h(i+1)/3.0;
    A(i+1,i+1) = A(i+1,i+1) - 1.0/h(i+1) + h(i+1)/3.0;
    A(i  ,i+1) = A(i  ,i+1) + 1.0/h(i+1) + h(i+1)/6.0;
    A(i+1,i  ) = A(i+1,i  ) + 1.0/h(i+1) + h(i+1)/6.0;
end;

% Create the load vector

b = zeros(N,1);

for i=1:N
    b(i,1) = (h(i)+h(i+1))*g(i+1)/3.0 + (h(i)*g(i)+h(i+1)*g(i+2))/6.0;
end;

% Compute the coefficients

c = A \ b;

% Compute the approximative solution on the plot-mesh

yp = interp1(t, [0; c; 0]', tp, 'linear');

maxer = max(abs(yp + tp - sin(tp) / sin(1.0)))