function [tp,yp] = leastsquares1(N,Np);
%
%   function [tp,yp] = leastsquares1(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 (1-t)*t^k.
%

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

% Create the coefficient matrix

A = zeros(N,N);

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

for ell=2:N
    k  = 1;
    kp = 2.0;
    lm = ell*(ell-1);
    lp = (ell+1)*ell;
    A(k,ell) = -kp*lm/(k+ell-2) + kp*lp/(k+ell-1) - kp/(k+ell) + kp/(k+ell+1) ...
               +lm/(k+ell-1) - lp/(k+ell) + 1.0/(k+ell+1) - 1.0/(k+ell+2) ...
               -lm/(k+ell) + lp/(k+ell+1) - 1.0/(k+ell+2) + 1.0/(k+ell+3);
    A(ell,k) = A(k,ell);
end;

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

% Create the right-hand side b

b = zeros(N,1);

for ell=1:N
    b(ell,1) = 1.0 + 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)))