function [tp,yp] = fem3(N,Np);
%
%   function [tp,yp] = fem3(N,Np);
%
% This function implements the finite element method with piecewise
% linear continuous basis functions for the PDE y'' + y = g(t), where
% g is the exponential function from page 3.90, 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.
%
% Here, the grid consists of N equidistant points from 0 to 0.1, and the
% two points 0.2 and 1.0.
%

t  = [linspace(0,0.1,N), 0.2, 1.0];
tp = linspace(0,1,Np);
h  = diff(t);

% Right-hand side evaluated at the mesh points

c1 = 100.0;
c2 = c1*c1;
g  = (4.0*c1*t - 2.0 - 2.0*c1) .* exp(-c1*t);

% Create the stiffness matrix

A = zeros(N,N);

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

for i=1:N-1
    A(i  ,i  ) = A(i  ,i  ) - 1.0/h(i+1) - c2 * h(i+1)/3.0;
    A(i+1,i+1) = A(i+1,i+1) - 1.0/h(i+1) - c2 * h(i+1)/3.0;
    A(i  ,i+1) = A(i  ,i+1) + 1.0/h(i+1) - c2 * h(i+1)/6.0;
    A(i+1,i  ) = A(i+1,i  ) + 1.0/h(i+1) - c2 * 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 .* (1.0-tp) .* exp(-c1*tp)))