function tumor_pde ( )

%*****************************************************************************80
%
%% tumor_pde solves the Chaplain and Stuart tumor angiogenesis PDE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2021
%
%  Author:
%
%    Original version by Chiara Villa.
%    Modifications by John Burkardt
%
%  Reference:
%
%    Mark Chaplain, Andrew Stuart,
%    A model mechanism for the chemotactic response of endothelial cells
%    to tumor angiogenesis factor,
%    IMA Journal of Mathematics Applied to Medicine & Biology,
%    Volume 10, Number 3, pages 149-168, September 1993.
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent advection-diffusion-reaction equations,
%    Springer, 2003,
%    ISBN: 978-3-662-09017-6
%
%    Chiara Villa,
%    Mathematical modelling of tumour-induced angiogensis,
%    University of St Andrews, 2018
%
  clc;
  clear ( 'all' );
  close ( 'all' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'tumor_pde:\n' );
  fprintf ( 1, '  MATLAB/Octave version "%s"\n', version ( ) );
  fprintf ( 1, '  Solve a tumor angiogenesis PDE.\n' );

  nx = 201;
  xmin = 0.0;
  xmax = 1.0;
  x = linspace ( xmin, xmax, nx );

  dt = 0.05;
  tmin = 0.0;
  tmax = 0.7;
  t = linspace ( tmin, tmax, 15 );
%
%  Solution requested at TSOLVE values.
%
  tsolve = [ 0.0, 0.1, 0.3, 0.5, 0.7 ];
%
%  M = 0 indicates "slab symmetry" of the region.
%  That is, 1D Cartesian coordinates with no symmetry imposed.
%
  m = 0;
%
%  PDEPE computes the solution.
%
  sol = pdepe ( m, @pdefun, @icfun, @bcfun, x, tsolve );
%
%  Display component 1 over time.
%
  figure ( ); 
  clf ( );
  hold ( 'on' );
    a1 = plot ( x, sol(1,:,1), 'linewidth', 3 );
    r1 = 't = 0.0';
    a2 = plot ( x, sol(2,:,1), 'linewidth', 3 );
    r2 = 't = 0.1';
    a3 = plot ( x, sol(3,:,1), 'linewidth', 3 );
    r3 = 't = 0.3';
    a4 = plot ( x, sol(4,:,1), 'linewidth', 3 );
    r4 = 't = 0.5';
    a5 = plot ( x, sol(5,:,1), 'linewidth', 3 );
    r5 = 't = 0.7';
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( 'Distance' );
  ylabel ( 'TAF concentration near tumor' );
  title ( 'Tumor Angiogenesis Factor' );
  legend ( [ a1; a2; a3; a4; a5 ], r1, r2, r3, r4, r5 );
  filename = 'taf.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Display component 2 over time.
%
  figure ( );
  clf ( );
  hold ( 'on' );
    b1 = plot ( x, sol(1,:,2), 'linewidth', 3 );
    s1 = 't = 0.0';
    b2 = plot ( x, sol(2,:,2), 'linewidth', 3 );
    s2 = 't = 0.1';
    b3 = plot ( x, sol(3,:,2), 'linewidth', 3 );
    s3 = 't = 0.3';
    b4 = plot ( x, sol(4,:,2), 'linewidth', 3 );
    s4 = 't = 0.5';
    b5 = plot ( x, sol(5,:,2), 'linewidth', 3 );
    s5 = 't = 0.7';
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( 'Distance' );
  ylabel ( 'Density' );
  title ( 'Endothelial cell density near tumor' );
  legend ( [ b1; b2; b3; b4; b5 ], s1, s2, s3, s4, s5, 'Location', 'North' );
  filename = 'density.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'tumor_pde:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

function [ c, f, s ] = pdefun ( x, t, u, DuDx )

%*****************************************************************************80
%
%% pdefun defines the PDEs.
%
%  Discussion:
%
%    pdepe() assumes we are treating one or more parabolic PDE's of the form
%
%      c(x,t,u,dudx) * du/dt = x^(-m) d/dx x^m f(x,t,u,dudx) + s(x,t,u,dudx)
%    
%    where for our problem, m = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2021
%
%  Author:
%
%    Original version by Chiara Villa.
%    Modifications by John Burkardt
%
%  Input:
%
%    real X, the spatial coordinate.
%
%    real T, the time coordinate.
%
%    real U(2), the solution at (X,T).
%
%    real DUDX(2), the spatial derivative of the solution at (X,T).
%
%  Output:
%
%    real C(2), the coefficients of dUdT
%
%    real F(2), the flux terms.
%
%    real S(2), the source terms.
%

%
%  Set physical parameters:
%
  alpha = 10.0;
  gamma = 1.0;
  lambda = 1.0;
  delta = 1.0;
  epsilon = 0.001;
  k = 0.75;
  mu = 100.0;
  beta = 4.0;
  cstar = 0.2;

  c = [ 1.0; 1.0 ];

  f = [ delta * DuDx(1); ...
       epsilon * DuDx(2) - k * u(2) * DuDx(1) ];

  g = max ( u(1) - cstar, 0.0 );

  s = [ - alpha * u(1) * u(2) / ( gamma + u(1) ) - lambda * u(2); ...
        mu * u(2) * ( 1.0 - u(2) ) * g - beta * u(2) ];

  return
end

function u0 = icfun ( x )

%*****************************************************************************80
%
%% icfun sets the initial condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2021
%
%  Author:
%
%    Original version by Chiara Villa.
%    Modifications by John Burkardt
%
%  Input:
%
%    real X, the spatial coordinate.
%
%  Output:
%
%    real U0(2), the initial condition for the solution at X.
%
  c0 = cos ( 0.5 * pi * x );
  if ( x < 1.0 )
    n0 = 0.0;
  else
    n0 = 1.0;
  end

  u0 = [ c0; n0 ];

  return
end

function [ pl, ql, pr, qr ] = bcfun ( xl, ul, xr, ur, t )

%*****************************************************************************80
%
%% bcfun sets the left and right boundary conditions.
%
%  Discussion:
%
%    The boundary conditions have the form
%
%      p(x,t,u) + q(x,t) * f(x,t,u,dudx) = 0
%
%  Input:
%
%    real XL, UL, the left spatial coordinate and solution value.
%
%    real XR, UR, the right spatial coordinate and solution value.
%
%    real T, the time coordinate.
%
%  Output:
%
%    real PL, QL, the coefficients for the left boundary condition.
%
%    real PR, QR, the coefficients for the right boundary condition.
%

%
%  On the left, U(1) = 1, and U(2) has zero flux.
%
  pl = [ ul(1)-1.0; 0.0 ];
  ql = [ 0.0; 1.0 ];
%
%  On the right, U(1) = 0, and U(2) = 1.
%
  pr = [ ur(1); ur(2)-1.0 ];
  qr = [ 0.0; 0.0 ];

  return
end

