%-----------------------------------------------------------------------
% Solve the diffusion-reaction problem via the finite difference approach.
%-----------------------------------------------------------------------
% The following is a one-dimensional diffusion-reaction equation.
%     Diffusion = Reaction
%
%   (d^2y/dx^2) = f(y)    BC: at x=0, dy/dx=0
%                             at x=1, y=1
%
% We devide x=[0,1] into n intervals.
% We approximate dy/dx     with     dy/dx = (y(i+1)-y(i-1))/(2*dx)
% We approximate d^2y/dx^2 with d^2y/dx^2 = (y(i+1)-2*y(i)+y(i-1))/dx^2
%
% The finite difference representation of the diffusion-reaction equation is:
%
%   (y(i+1)-2*y(i)+y(i-1))/dx^2 - f(y(i)) = 0         for  i=2,3,...,n
%   g(i) = y(i+1)-2*y(i)+y(i-1) - dx*dxf(y(i)) = 0    for  i=2,3,...,n
%
% We solve g(y)=0 with, say, the Newton's method.
%
% An example with n intervals
% i=1:  g(1) = y(1)-y(2) = 0
% i=2:  g(2) = y(3)-2*y(2)+y(1) - dx*dx*f(y(2)) = 0
% i=3:  g(3) = y(4)-2*y(3)+y(2) - dx*dx*f(y(3)) = 0
% i=4:  g(4) = y(5)-2*y(4)+y(3) - dx*dx*f(y(4)) = 0
%  :                                                                         :
% i=n:  g(n) = y(n+1)-2*y(n)+y(n-1) - dx*dx*f(y(n)) = 0
% i=n+1: g(n+1) = y(n+1) - 1 = 0
%
%-----------------------------------------------------------------------
% 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(' ')

% Divide x=[0,1] into n subintervals -----------------------------------
      n = 100;      % n=number of subintervals
      dx = 1/n;
      for i=2: n+1
        x(i) =  dx*(i-1);
      end

% Function for saturation (Michaelis-Menten) kinetics
%     f     = alpha*y(i)/(beta+y(i))
      alpha = 1;
      beta  = 1;

% Provide initial guess ------------------------------------------------
      y = ones(n+1,1);

% Iterate with Newton's Method -----------------------------------------
    for k=1: 3
      g = zeros(n+1,1);
      dgdy = zeros(n+1,n+1);

% Assign vector g and Jacobian matrix dgdy -----------------------------
      for i=2: n
        f    = alpha*y(i)/(beta+y(i));
        dfdy = alpha*beta/(beta+y(i))^2;
        g(i) = y(i+1) - 2*y(i) + y(i-1) - dx*dx*f;
        dgdy(i,i-1) = 1;
        dgdy(i,i)   = -2 -dx*dx*dfdy;
        dgdy(i,i+1) = 1;
      end

% B.C. -----------------------------------------------------------------
% At x=0, dy/dx=0
      g(1)       = y(1)-y(2);
      dgdy(1,1)  = 1;
      dgdy(1,2)  = -1;
% At x=1, y=1
      g(n+1)     = y(n+1)-1;
      dgdy(n+1,n+1) = 1;

% Solve the linear set of equations with predivision -------------------
      y = y-dgdy\g;
% Equally valid ways of finding solution to a linear system of equations.
%     y = y-inv(A)*b;
%     y = y-A^(-1)*b;
    end

% 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))])