function [tp,yp] = leastsquares2(N,Np);
%
%   function [tp,yp] = leastsquares2(N,Np);
%
% This function implements the least squares 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 sin(k*pi*t).
%

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

% Create the coefficient matrix

A = zeros(N,N);

for k=1:N
    A(k,k) = (1.0 - k*k*pi*pi)*(1.0 - k*k*pi*pi)/2.0;
end;

% Create the right-hand side b

b = zeros(N,1);

for ell=1:N
    b(ell,1) = (-1)^ell * (1.0/(ell*pi) - ell*pi);
end;

% Solve for the coefficients

c = A \ b;

% Compute the function values

yp = zeros(Np,1);

for k=1:N
    yp = yp + c(k) * sin(k*pi*tp);
end;

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