%  Math 581, Fall 2006.
%  MATLAB code written by Curt Vogel.
%
%  Compute finite difference approximation to solution of ODE BVP
%      -d^2u/dx^2 = f(x), 0 < x < 1,
%       u(0) = 1,
%       u(1) = 0.
%  Test the code by taking f(x) = 9*pi^2/4 * cos(3*pi*x/2), with  
%  exact solution u_exact = cos(3*pi*x/2).
%
%    Approximate solutions are computed on a sequence of grids.
%  For each grid the solution error is computed. A log-log
%  plot of solution error vs. grid spacing is generated.

  nlevels = input(' No. of grid levels = ');
  nvec = zeros(nlevels,1);
  for k = 1:nlevels
    N = 2 * 2^k - 1        %  No. of unknowns = no. of interior grid points.
    nvec(k) = N;           %  Save no. of unknowns.
    h = 1/(N+1)            %  Grid spacing.
    x = [0:h:1]';          %  Grid points, including boundary.
    x_interior = x(2:N+1); %  Interior grid points.
  
%  Set up N-by-N matrix A_h = tridiag(-1,2,-1). Use MATLAB's sparse
%  storage mode.  First set up main diagonal and subdiagonal of A_h.
%  Since A_h is symmetric, the sub- and superdiagonals are the same.
  
    Ah_diag = 2 * ones(N,1);
    Ah_sub = -1 * ones(N-1,1);
    A_h = spdiags([[Ah_sub;0],Ah_diag,[0;Ah_sub]],[-1 0 1],N,N);

%  Set up right-hand-side vector b_h. First evaluate r.h.s function f(x)
%  at interior grid points. Then add contribution from boundary conditions.

    rhs_function = 9*pi^2/4 * cos(1.5*pi*x_interior);
    u_left = 1;   %  Left BC
    u_right = 0;  %  Right BC
    b_h = h^2 * rhs_function;
    b_h(1) = b_h(1) + u_left;
    b_h(N) = b_h(N) + u_right;
    
%  Solve linear system A_h*uh_interior = b_h and append boundary values. 
    
    uh_interior = A_h\b_h;
    uh = [u_left; uh_interior; u_right];
    
%  Evaluate exact solution and compute discrete solution error.

    u_exact = cos(3*pi*x/2);
    Eh = norm(uh - u_exact,'inf');
    
%  Save error and mesh spacing.
    
    Ehvec(k) = Eh;
    hvec(k) = h;

     (u_exact.-uh(ii))./uh

%  Plot true and approximate solutions in figure 1.

  figure(1)
      plot(x,u_exact,'*-', x,uh,'o-')
      xlabel('x axis')
      ylabel('u(x)')
      title_str = ['Grid spacing h = ', num2str(h)];
      title(title_str)
     # legend('True Solution','FD Approximation')

      disp(' Hit any key to continue.'); 
      pause

  end

%  Plot error vector Eh vs. h on a log-log scale.

  figure(2)
    loglog(hvec,Ehvec,'o', hvec,Ehvec)
    xlabel('mesh size h')
    ylabel('E_h = max |u_h - u_{exact}|')
    title('Log-log plot of solution error E_h vs. mesh size h')
        
%  Estimate rate of convergence.

  conv_rate = mean( diff(log(Ehvec)) ./ diff(log(hvec)) );
  fprintf('Estimated rate of convergence is h^p, where p = %4.3e\n',conv_rate);