function diffusion_pde_ode45 ( )

%*****************************************************************************80
%
%% diffusion_pde_ode45() solves the diffusion PDE using ode45().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 September 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent diffusion-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'diffusion_pde_ode45():\n' );
  fprintf ( 1, '  Solve the diffusion PDE in 1D,\n' );
  fprintf ( 1, '    du/dt - mu d2u/dx2 = 0\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with zero Neumann boundary conditions:\n' );
  fprintf ( 1, '    u''(0) = 0\n' );
  fprintf ( 1, '    u''(1) = 0\n' );
  fprintf ( 1, '  and diffusion coefficient\n' );
  fprintf ( 1, '    mu = constant\n' );
  fprintf ( 1, '  and initial condition\n' );
  fprintf ( 1, '    u(0,x) = (10x-6)^2 (8-10x)^2 for 0.6 <= x <= 0.8\n' );
  fprintf ( 1, '           = 0 elsewhere.\n' );
  fprintf ( 1, '  and NX equally spaced nodes in X,\n' );
  fprintf ( 1, '  and NT equally spaced points in T,\n' );
  fprintf ( 1, '  using ode45().\n' );

  [ mu, t0, tstop, xmin, xmax ] = diffusion_parameters ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Parameters:\n' );
  fprintf ( 1, '    mu    = %g\n', mu );
  fprintf ( 1, '    t0    = %g\n', t0 );
  fprintf ( 1, '    tstop = %g\n', tstop );
  fprintf ( 1, '    xmin  = %g\n', xmin );
  fprintf ( 1, '    xmax  = %g\n', xmax );

  nx = 101;
  x = linspace ( xmin, xmax, nx );

  u0 = diffusion_initial ( x );
%
%  Call ode45() to compute the solution.
%
  f = @(t,u)diffusion_deriv(t,u,x);
  tspan = [ t0, tstop ];

  [ t, u ] = ode45 ( f, tspan, u0 );

  nt = length ( t );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of adaptive time steps = %d\n', nt );
%
%  Initialize the movie.
%
  colormap ( 'jet' );
  moviename = 'diffusion_pde_ode45.avi';
  Video = VideoWriter ( moviename, 'Motion JPEG AVI' );
  Video.FrameRate = 10;
  open ( Video );
%
%  Compute the solution at each time.
%
  umax = max ( max ( u ) );
  umin = min ( min ( u ) );

  for j = 1 : nt
%
%  Display the solution.
%
    clf ( );
    hf = figure ( 1 );

    plot ( x, umin, 'r-', ...
           x, umax, 'r-', ...
           x, u(j,:), 'b-', ...
           'linewidth', 3 );

    grid ( 'on' );
    xlabel ( '<--- x --->' );
    ylabel ( '<--- u(x,t) --->' );
    axis ( [ 0.0, 1.0, -0.05, 1.05 ] )
    title_string = sprintf ( 'diffusion PDE ode45, Step %d, Time = %8.4f\n', j, t(j) );
    title ( title_string );

    if ( j == 1 )
      filename = 'diffusion_pde_ode45_initial.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( j == nt )
      filename = 'diffusion_pde_ode45_final.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    end

    frame = getframe ( 1 );
    writeVideo ( Video, frame );

  end
%
%  Close the movie.
%
  close ( Video );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );
  return
end

