%
% bvp_2c.m 
%
% second order finite difference method for the general ODE bvp
%   a(x)u''(x)+b(x)u'(x)+c(x)u(x) = f(x), 
%   u(ax)=alpha, u(bx)=beta 
%
% From  http://www.amath.washington.edu/~rjl/fdmbook/  (2007)

ax = 0;
bx = 1;
alpha = 1;    
beta = exp(1);     % Dirichlet boundary condtion at bx

utrue = @(x) exp(x);  % true soln

afun=@(x) 1+x;
bfun=@(x) x;
cfun=@(x) exp(x);
%construct f(x) based on the ODE
f = @(x) afun(x).*exp(x)+bfun(x).*exp(x)+cfun(x).*exp(x);  % right hand side function


% true solution on fine grid for plotting:
xfine = linspace(ax,bx,101);
ufine = utrue(xfine);

% Solve the problem for ntest different grid sizes to test convergence:
m1vals = [10 20 40 80 160 320];
ntest = length(m1vals);
hvals = zeros(ntest,1);  % to hold h values
E = zeros(ntest,1);      % to hold errors

for jtest=1:ntest
  m1 = m1vals(jtest);
  m2 = m1 + 1;
  m = m1 - 1;                 % number of interior grid points
  hvals(jtest) = (bx-ax)/m1;  % average grid spacing, for convergence tests

  % set grid points:  
  gridchoice = 'uniform';          % see xgrid.m for other choices
  x = xgrid(ax,bx,m,gridchoice);   

  % set up matrix A (using sparse matrix storage):
  A = spalloc(m2,m2,3*m2);   % initialize to zero matrix

  % first row for Neumann BC at ax:
  A(1,1:3) = fdcoeffF(0, x(1), x(1:3)); 

  % interior rows:
  % a(x)u''(x)+b(x)u'(x)+c(x)u(x)=f(x)
  for i=2:m1
     A(i,i-1:i+1) = afun(x(i))*fdcoeffF(2, x(i), x((i-1):(i+1)))... %a(x)u''(x)
         +bfun(x(i))*fdcoeffF(1, x(i), x((i-1):(i+1)))... %b(x)u'(x)
         +[0 cfun(x(i)) 0]; %c(x)u(x), last term only multiply u(x(i))
  end

  % last row for Dirichlet BC at bx:
  A(m2,m:m2) = fdcoeffF(0,x(m2),x(m:m2)); 
  
  % Right hand side:
  F = f(x); 
  h=hvals(jtest);
  F(1) = alpha;   
  F(m2) = beta;
  
  % solve linear system:
  U = A\F;


  % compute error at grid points:
  uhat = utrue(x);
  err = U - uhat;
  E(jtest) = max(abs(err));  
  disp(' ')
  disp(sprintf('Error with %i points is %9.5e',m2,E(jtest)))

  clf
  plot(x,U,'o')  % plot computed solution
  title(sprintf('Computed solution with %i grid points',m2));
  hold on
  plot(xfine,ufine)  % plot true solution
  hold off
  
  % pause to see this plot:  
  %drawnow
  %input('Hit <return> to continue ');
  
  end

error_table(hvals, E);   % print tables of errors and ratios
error_loglog(hvals, E);  % produce log-log plot of errors and least squares fit

