function [ t, y ] = leapfrog ( dydt, tspan, y0, n )

%*****************************************************************************80
%
%% leapfrog() solves an ODE using the leapfrog method.
%
%  Discussion:
%
%    There are two crucial assumptions:
%
%    1) The number of variables is 2.
%    2) Written as a first order system, the ODE has the form
%       u' = v
%       v' = f(t,u)   NO DEPENDENCE ON V
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 April 2020
%
%  Author:
%
%    Original Python version by John D Cook;
%    This version by John Burkardt
%
%  Reference:
%
%    John D Cook,
%    Leapfrog integrator,
%    https://www.johndcook.com/blog/2020/07/13/leapfrog-integrator/
%    13 July 2020.
%
%  Input:
%
%    function dydt: evaluates the right hand side of the ODE.
%
%    real tspan(2): contains the initial and final times.
%
%    real y0(2): the initial condition.
%
%    integer n: the number of steps to take.
%
%  Output:
%
%    real t(n+1), y(n+1,2): the times and solution values.
%
  t0 = tspan(1);
  tstop = tspan(2);
  dt = ( tstop - t0 ) / n;

  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, 2 );

  for i = 1 : n + 1

    if ( i == 1 )
      t(1)   = t0;
      y(1,1) = y0(1);
      y(1,2) = y0(2);
      anew   = dydt ( t, y(i,:) );
    else
      t(i)   = t(i-1) + dt;
      aold   = anew;
      y(i,1) = y(i-1,1) + dt * ( y(i-1,2) + 0.5 * dt * aold(2) );
      anew   = dydt ( t, y(i,:) );
      y(i,2) = y(i-1,2) + 0.5 * dt * ( aold(2) + anew(2) );
    end

  end

  return
end

