% P1 
% Checks out bandLU and BandSolve

  randn('seed',0)

% Small example for correctness check...
  n = 5;
  a = randn(n-1,1);
  b = randn(n,1);
  [c,d,f] = BandLU(a);
  x = BandSolve(c,d,f,b);
  clc
  disp(['x =' sprintf('%10.6f ',x)])
  disp(['c =' sprintf('%10.6f ',c)])
  disp(['d =' sprintf('%10.6f ',d)])
  disp(['f =' sprintf('%10.6f ',f)])
 
% Let's solve  u + u' = 10exp(-x)cos(10x) across [0,2*pi]
% Solution: u(x) = exp(-x)*sin(10x)

  n = 1000;
  x = linspace(0,pi,n)';
  h = pi/(n-1);
% The discrete equation at x(k) is  
%            u(k) + (u(k+1)-u(k-1)/2h = exp(-x(k))*sin(10x(k))
  a = (1/(2*h))*ones(n-3,1);
  b = 10*exp(-x(2:n-1)).*cos(10*x(2:n-1));
  [c,d,f] = BandLU(a);
  u = BandSolve(c,d,f,b);
  plot(x,[0;u;0],x,0*x,'--r')
  axis([0 pi -1 1])
  shg