function [ u, x, t ] = kdv_mid ( nx, nt )
%
% function [ u, x, t ] = kdv_mid ( nx, nt )
%
%
%  Compute a solution of the KdV wave equation, using symmetric approximations
%  to the spatial derivatives, and the midpoint method for time integration.
%
%  Discussion:
%
%    The spatial domain is assumed to extend from -6 to 6.  Increasing NX
%    decreases the spatial stepsize but does not increase the size of the domain.
%
%    The time computation is done in fixed steps of size 0.002.  Increasing
%    NT does not change the stepsize, but increases the final time at which
%    results are computed.
%
%  Modified:
%
%    10 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NT, the number of points in the X direction, and
%    the number of time values in the time direction.  Reasonable results
%    can be obtained with NX = 65 and NT = 121.
%
%    Output, real U(NT,NX), the computed solution.
%
%    Output, real X(NX), T(NT), the gridpoint values of space and time.
%
xmin = -6.0;
xmax =  6.0;
dx = ( xmax - xmin ) / ( nx - 1 );
for ix = 1 : nx
  x(ix) = ( ( nx - ix ) * xmin + ( ix - 1 ) * xmax ) / ( nx - 1 );
end

dt = 0.002;
for it = 1 : nt
  t(it) = ( it - 1 ) * dt;
end

u = zeros (nt,nx);
%
%  Set the initial condition.
%  N_SOLITON is the number of solitons we want to see.  You can increase
%  the value to 3 or 4, but the code isn't ready for significantly higher values.
%
N_SOLITON = 2;

it = 1;
for ix = 1 : nx
  u(it,ix) = N_SOLITON * ( N_SOLITON + 1 ) * ( sech ( x(ix) ) )^2;
end
%
%  First step uses Euler's method.
%
it = 2;

  u(it,1) = u(it-1,1);
  u(it,2) = u(it-1,2);

  for ix = 3 : nx-2
    u(it,ix) = u(it-1,ix) + dt * (...
      - 6.0 * u(it-1,ix) * ( u(it-1,ix+1) - u(it-1,ix-1) ) / ( 2.0 * dx ) ...
      - ( u(it-1,ix+2) - 2 * u(it-1,ix+1) + 2 * u(it-1,ix-1) - u(it-1,ix-2) ) / ( 2 *dx^3 ) );
  end

  u(it,nx-1) = u(it-1,nx-1);
  u(it,nx) = u(it-1,nx);
%
%  Subsequent steps use the midpoint method.
%
for it = 3 : nt

  u(it,1) = u(it-1,1);
  u(it,2) = u(it-1,2);

  for ix = 3 : nx-2
    u(it,ix) = u(it-2,ix) + 2.0 * dt * (...
      - 6.0 * u(it-1,ix) * ( u(it-1,ix+1) - u(it-1,ix-1) ) / ( 2.0 * dx ) ...
      - ( u(it-1,ix+2) - 2 * u(it-1,ix+1) + 2 * u(it-1,ix-1) - u(it-1,ix-2) ) / ( 2 *dx^3 ) );
  end

  u(it,nx-1) = u(it-1,nx-1);
  u(it,nx) = u(it-1,nx);

end