%
% chap2page35.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
% Using central differences on a uniform grid.
%
% Solve ODE (2.64-2.65) with FM system (2.68) 
% Based on  http://www.amath.washington.edu/~rjl/fdmbook/chapter2  (2007)
set(0, 'defaultaxesfontsize',16,'defaultaxeslinewidth',1,...
    'defaultlinelinewidth',1,'defaultpatchlinewidth',0.7,...
    'defaulttextfontsize',16);

%setup the ODE 
ax = 0;
bx = 1;
alpha = 1;
beta = exp(1);
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


% 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;
  h = (bx-ax)/m1;
  hvals(jtest) = h;
  x = linspace(ax,bx,m2)';   % m+2 points counting the boundaries

  % set up matrix A (using sparse matrix storage):
  % construct system (2.67) by first compute 3 main diagonals
  %   
  supdiag=afun(x)+h*bfun(x)/2;
  maindiag=h^2*cfun(x)-2*afun(x);
  subdiag=afun(x)-h*bfun(x)/2;
  
  %In spdiags syntax, if a column of B is longer than the diagonal it is replacing, 
  %spdiags takes elements of super-diagonals from the lower part of the column of B, 
  %and elements of sub-diagonals from the upper part of the column of B.
  %test this code: e=[1:5]';full(spdiags([e e e], [-1 0 1], 5, 5))
  B = [subdiag(3:end) maindiag(2:end-1) supdiag(1:end-2)]; %pay attention to shifted index
  A = 1/h^2 * spdiags(B, [-1 0 1], m, m);
  
  % Right hand side: see (2.68)
  F = f(x(2:end-1)); %only interial points
  F(1) = F(1)-(afun(x(2))/h^2-bfun(x(2))/2/h)*alpha ;  % Approach 1 or 2 from Section 2.12
  F(m) = F(m)-(afun(x(end-1))/h^2+bfun(x(end-1))/2/h)*beta;
  
  % solve linear system:
  U = A\F;

  % compute errors:
  uhat = utrue(x(2:end-1));
  err = U - uhat;
  E(jtest) = max(abs(err)); 
 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

figure
%plot the solutions and the error curves
subplot(2,1,1); hold on
plot(x,[alpha;U;beta],'--')
plot(x,utrue(x),'-r');
legend('FD solution','True solution');
xlabel('x'); ylabel('{u_h(x)} and {u^*(x)}');
title('solution curves')

subplot(2,1,2); hold on
semilogy(x,[0;err;0])
title('error curve')
xlabel('x'); ylabel('{u_h(x)-u^*(x)}');
