% solution for HW4, #3, especially part c) (see end).  Barnett 10/28/08
% Makes use of gauss.m code.

% define RHS and kernel functions. interval [0,1]
f = @(t) exp(t) + (exp(t+1)-1)./(t+1);    % important for it to handle vectors
k = @(s,t) exp(s.*t);
uex = @(t) exp(t);                        % exact soln

Ns = 2:10;                                % for Gauss in a) ii)
%Ns = 10:10:300;                          % for comp trap in a) i)
for i=1:numel(Ns)
  N = Ns(i);
  [x w] = gauss(N);  x = (x+1)/2; w = w/2; % gauss, rescaled to interval [0,1]
  %x = (0:N-1)'/(N-1); w = [.5 ones(1,N-2) .5]/(N-1); % comp trap, likewise
  b = f(x);                               % RHS col vec
  [ss tt] = meshgrid(x);
  A = repmat(w, [numel(x) 1]) .* k(ss,tt); % multiply each col by resp. weight
  u = (eye(numel(x)) + A)\b;
  fprintf('N = %d: cond=%g, error sup norm = %g\n', N, ...
          cond(eye(numel(x))+A), max(abs(u - uex(x))))  % note I didn't plot
end

% c) plot the answer using the Nystrom interpolation formula:
N = 5; [x w] = gauss(N);  x = (x+1)/2; w = w/2;    % Gauss for N=5, rescaled
b = f(x);
[ss tt] = meshgrid(x);
A = repmat(w, [numel(x) 1]) .* k(ss,tt);
uj = (eye(numel(x)) + A)\b;                % as above, soln vec, u^(N) at nodes
t = 0:.01:1;                               % grid to plot u on, NOT just nodes!
[ss tt] = meshgrid(x, t);                  % notice rectangular (tall)
At = repmat(w, [numel(t) 1]) .* k(ss,tt);  % eval matrix: mult cols by weights
u =  -At*uj + f(t)';                       % the (N) from lect, w/ RHS col vec
figure; plot(t, u' - uex(t), '-'); hold on; plot(x, 0, 'k.');
xlabel x;ylabel('u^{(N)} - u'); title('error function for Nystrom solution');