function pendulum_double_ode45 ( )

%*****************************************************************************80
%
%% pendulum_double_ode45() uses ode45() on pendulum_double_ode().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 May 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'pendulum_double_ode45():\n' );
  fprintf ( 1, '  Use ode45() on pendulum_double_ode().\n' );
%
%  Get the parameters.
%
  [ g, m1, m2, l1, l2, t0, y0, tstop ] = pendulum_double_parameters ( );
%
%  Reset the default value of TSTOP.
%
  tstop = 50.0;
  [ g, m1, m2, l1, l2, t0, y0, tstop ] = pendulum_double_parameters ( ...
    g, m1, m2, l1, l2, t0, y0, tstop );
%
%  Report parameter values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  parameters:\n' );
  fprintf ( 1, '    g     = %g\n', g );
  fprintf ( 1, '    m1    = %g\n', m1 );
  fprintf ( 1, '    m2    = %g\n', m2 );
  fprintf ( 1, '    l1    = %g\n', l1 );
  fprintf ( 1, '    l2    = %g\n', l2 );
  fprintf ( 1, '    t0    = %g\n', t0 );
  fprintf ( 1, '    y0    = (%g,%g,%g,%g)\n', y0 );
  fprintf ( 1, '    tstop = %g\n', tstop );

  f = @ pendulum_double_deriv;
  tspan = [ t0, tstop ];

  [ t, theta ] = ode45 ( f, tspan, y0 );

  n = size ( t, 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of variable sized steps = %d:\n', n );

  x1 =   l1 * sin ( theta(:,1) );
  y1 = - l1 * cos ( theta(:,1) );

  x2 = x1 + l2 * sin ( theta(:,3) );
  y2 = y1 - l2 * cos ( theta(:,3) );
%
%  Solution plot.
%
  figure ( 1 );
  clf ( );
  plot ( t, x1, 'g', ...
         t, y1, 'r', ...
         'LineWidth', 2 );
  grid ( 'on' );
  xlabel ( 'Time' );
  ylabel ( 'x1(t), y1(t)' );
  legend ( 'x1(t)', 'y1(t)' );
  title ( 'pendulum_double_ode(): ode45(): (x1,y1)' )
  filename = 'pendulum_double_ode45_x1y1.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  figure ( 2 );
  clf ( );
  plot ( t, x2, 'g', ...
         t, y2, 'r', ...
         'LineWidth', 2 );
  grid ( 'on' );
  xlabel ( 'Time' );
  ylabel ( 'x2(t), y2(t)' );
  legend ( 'x2(t)', 'y2(t)' );
  title ( 'pendulum_double_ode(): ode45(): (x2,y2)' )
  filename = 'pendulum_double_ode45_x2y2.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Phase plane plot.
%
  figure ( 3 );
  clf ( );
  hold ( 'on' );
  plot ( x1, y1, 'r-', 'linewidth', 3 );
  plot ( x2, y2, 'b-', 'linewidth', 1 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- x(t) --->' );
  ylabel ( '<--- y(t) --->' );
  legend ( '(x1,y1)', '(x2,y2)' );
  axis ( 'equal' );
  title ( 'pendulum_double_ode(): ode45(): phase' )
  filename = 'pendulum_double_ode45_phase.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Energy plot.
%
  h = pendulum_double_conserved ( theta );

  figure ( 4 );
  clf ( )
  hold ( 'on' );
    plot ( t, h, 'r-', 'linewidth', 3 );
    plot ( tspan, [h(1),h(1)], 'k:', 'linewidth', 2 );
    plot ( tspan, [0.0,0.0], 'k-', 'linewidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- energy --->' );
  title ( 'pendulum_double_ode(): ode45(): conservation plot' )
  filename = 'pendulum_double_ode45_conservation.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

