function [tp,yp] = galerkin1(N,Np);
%
%   function [tp,yp] = galerkin1(N,Np);
%
% This function implements the Galerkin method
% for the problem y'' + y = -t subject to homogeneous
% Dirichlet boundary conditions. N denotes the number 
% of basis functions, and Np denotes the number of
% points in the output arrays for xp and up, for plotting
% the solution. The basis functions are (1-t)*t^k.
%

tp = linspace(0,1,Np)';

% Create the coefficient matrix

A = zeros(N,N);

k   = 1;
ell = 1;
kp  = (k+1)*k;
A(ell,k) = - kp/(k+ell) ...
           + 1.0/(k+ell+1) - 1.0/(k+ell+2) ...
           + kp/(k+ell+1) ...
           - 1.0/(k+ell+2) + 1.0/(k+ell+3);

for k=1:N
    for ell=1:N
        if (k*k+ell*ell>1)
            km = k*(k-1);
            kp = (k+1)*k;
            A(ell,k) =   km/(k+ell-1) - kp/(k+ell) ...
                       + 1.0/(k+ell+1) - 1.0/(k+ell+2) ...
                       - km/(k+ell) + kp/(k+ell+1) ...
                       - 1.0/(k+ell+2) + 1.0/(k+ell+3);
        end;
    end;
end;

% Create the right-hand side b

b = zeros(N,1);

for ell=1:N
    b(ell,1) = 1.0/(ell+3.0) - 1.0/(ell+2.0);
end;

% Solve for the coefficients

c = A \ b;

% Compute the function values

yp = c(N) * ones(Np,1);

for k=N-1:-1:1
    yp = yp .* tp + c(k);
end;

yp = yp .* tp .* (1.0 - tp);

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