%-----------------------------------------------------------------------
% Solve the diffusion-reaction problem via the finite difference approach.
%   (d^2y/dx^2) = f(y)    BC: at x=0, dy/dx=0
%                             at x=1, y=1
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Output a program header ----------------------------------------------
      disp('-------------------------------------------------------')
      disp('Find concentration profile of a diffusion-reaction eqn.')
      disp('Reaction rate = alpha*y/(beta+y)                       ')
      disp('-------------------------------------------------------')
      disp(' ')

% Find the concentration profile ---------------------------------------
      alpha=input('Enter alpha: ');
      beta =input('Enter beta:  ');
      [x, y] = finitedb(alpha, beta);
      n = max(size(x))-1   % number of intervals

% Output the answer (table) --------------------------------------------
      disp('The concentration profile is ')
      disp('  ------------------------------')
      disp('  Point#  Position Concentration')
      disp('  ------------------------------')
      for i=1:n+1
        fprintf('%6i%10.2f%12.4f\n', i, x(i), y(i) )
      end
      disp('  -----------------------------')

% Output the answer (graph) --------------------------------------------
      plot(x,y)
      xlabel('Position')
      ylabel('Concentration')
      title('Concentration Profile in Diffusion-Reaction System')
      text(0.1, 0.95, [ num2str(n) ' intervals'])
      text(0.05, y(1)+0.02, [ 'y(0)= ' num2str(y(1))])