function [ xp, yp ] = ode_midpoint_system_VS ( f_string, a, b, ya, tau0, tol )

%*****************************************************************************80
%
%% ode_midpoint_system_VS: variable step midpoint method for an ODE system.
%
%  Discussion:
%
%    This function uses the midpoint method to estimate the solution
%    to a system of ODE's y'=f(x,y) from a to b.  A variable step size
%    TAU is used, adjusted so that the local truncation error is 
%    about tolerance TOL, 
%
%    The user enters the name of an M file that evaluates f(x,y), the values 
%    of a and b, the initial condition vector for y at a, 
%    the initial time step TAU0 and tolerance TOL.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%      [ xp, yp ] = ode_midpoint_system_VS ( 'lotka(x,y)', 0.0, 250.0, [6.0,2.0], 25.e-3 , 1.e-4 );
%      [ xp, yp ] = ode_midpoint_system_VS ( 'lotka(x,y)', 0.0, 250.0, [2.0,3.0], 1.e-4 , 1.e-6 );
%      [ xp, yp ] = ode_midpoint_system_VS ( 'simplependulum(x,y)', 0.0, 200, [0.9*pi,0.0], 1.e-4 , 1.e-6 );
%    or the arguments can be prompted for:
%      [ xp, yp ] = ode_midpoint_system ( );
%      Enter function M-file name, such as lotka(x,y): lotka(x,y)
%      Enter starting point a: 0.0
%      Enter ending point b: 10.0
%      Enter the initial condition vector y(a), such as [1,2,3]: [4,2]
%      Enter the initial timesteps: 1.e-3
%      Enter tolerance for LTE: 1.e-5
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle f: the name of the M file that evaluates the
%    right hand side of the differential equation.  
%
%    real a, b: the starting and ending points.
%
%    real ya(m): the initial conditions at a. 
%
%    real tau0: the initial timestep to take.
%
%    real tol: the tolerance for the local truncation error.
%
%  Output:
%
%    real xp(n+1), yp(m,n+1): the sequence of solution estimates, 
%    which can be plotted.
%

%
%  Get f, a, b, ya, tau0, tol.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function M-file name, such as lotka(x,y): ', 's' );
  end
  f_element = vector_to_element ( f_string );
  f = str2func ( [ '@(x,y)', f_element ] );

  if ( nargin < 2 )
    a = input ( 'Enter the starting point a: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end

  if ( nargin < 3 )
    b = input ( 'Enter the ending point b: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end

  if ( nargin < 4 )
    ya = input ( 'Enter the initial condition vector y(a), such as [1,2,3]: ' );
  elseif ( ischar ( ya ) )
    ya = str2num ( ya );
  end
  ya = ya(:);
  m = length ( ya );

  if ( nargin < 5 )
    tau0 = input ( 'Enter the initial timestep tau0: ' );
  elseif ( ischar ( tau0 ) )
    tau0 = str2num ( tau0 );
  end

  if ( nargin < 6 )
    tol = input ( 'Enter the tolerance for LTE: ' );
  elseif ( ischar ( tol ) )
    tol = str2num ( tol );
   end
  
  theta = 0.500;

  hold all
  xp=[a] ;                         % t0
  yp(1,:) = ya(:) ;                % y0
  taunew = tau0;                  % tau0
  tau=[taunew];
  iter = 1;
%
% compute y1,  
%
  xm = xp(iter) + theta * taunew; % t1/2
  ym(1,:) = yp(iter,:);
    for j = 1 : 10
      ym(1,:) = yp(iter,:) + theta * taunew * ( f(xm,ym(1,:)) )';
    end
  yp(iter+1,:) = (       1.0 / theta ) * ym(1,:) ...
              + ( 1.0 - 1.0 / theta ) * yp(iter,:);       % y1 = yp(2,:)
  xp(iter+1)=xp(iter) + taunew ;                          % t1 =xp(2)
  tau(iter)=taunew;
  iter = iter+1;                                          % iter = 2 
%
% compute y2,  
%
  xm = xp(iter) + theta * taunew; % t3/2
  ym(1,:) = yp(iter,:);
    for j = 1 : 10
      ym(1,:) = yp(iter,:) + theta * taunew * ( f(xm,ym(1,:)) )';
    end
  yp(iter+1,:) = (       1.0 / theta ) * ym(1,:) ...
              + ( 1.0 - 1.0 / theta ) * yp(iter,:);        % y2 = yp(3,:)
  xp(iter+1)=xp(iter) + taunew;                            % t2 = xp(3)
  tau(iter)=taunew;
  iter = iter+1;                                           % iter = 3 
% so far tau0=tau1

tnew = xp(end);
taunew = tau(end);
count = 2; 

while tnew <=b
    count=count+1;
    tau(iter) = taunew;                                     % tau2 ~ tau(3)
    % evaluate y(n+1)
      xm = xp(iter) + theta * taunew; % tn+1/2
      ym(1,:) = yp(iter,:);
    for j = 1 : 10
      ym(1,:) = yp(iter,:) + theta * taunew * ( f(xm,ym(1,:)) )';
    end
      yp(iter+1,:) = (       1.0 / theta ) * ym(1,:) ...       % y3~y(4)
              + ( 1.0 - 1.0 / theta ) * yp(iter,:) ;          % yn = yp(n+1,:)
    % evaluate LTE using AB2-like  
        %first evaluate the AB2-like solution
        uAB2(iter,:) = yp(iter,:) * (tau(iter)+tau(iter-1))*(tau(iter)+tau(iter-1)+tau(iter-2)) / (tau(iter-1)*(tau(iter-1)+tau(iter-2)))...
            - yp(iter-1,:) * tau(iter)*(tau(iter)+tau(iter-1)+tau(iter-2)) / (tau(iter-1)*tau(iter-2))...
            + yp(iter-2,:) * tau(iter)*(tau(iter)+tau(iter-1)) / (tau(iter-2)*(tau(iter-1)+tau(iter-2)));
         % evaluate the (variable) error coefficient in the LTE R_n:
        R = 1/24+ 1/8*(1+tau(iter-1)/tau(iter))*(1+2*tau(iter-1)/tau(iter) + tau(iter-2)/tau(iter));
         % the approximation of the LTE ^T_{n+1}:
        Tnplus1(iter,:) = ( yp(iter+1,:) - uAB2(iter,:) )*1/(1-1/(24*R) );
        taunew = 0.95*tau(end)*abs(tol/norm(Tnplus1(iter,:)))^(1/3);
%        fprintf('The temporary timestep is %d\n',taunew);
        if norm(Tnplus1(iter,:))<= tol
%            fprintf('The norm of the LTE %f is less than the tolerance %f\n',norm(Tnplus1(iter,:)),tol);
            xp(iter+1) = xp(iter) + taunew ;                          % tn = xp(n+1)
            tnew=xp(end);
 %           fprintf('The current time is %d\n',tnew);
 %           fprintf('The new timestep is %d\n',taunew);      
            iter = iter+1;                                           % iter = ... 
        end
%      fprintf('The norm of the LTE %e is more than the tolerance %e\n',norm(Tnplus1(end,:)),tol);  
end
count
iter
fprintf('The number of time intervals is %d\n',iter(end));
fprintf('The minimum timestep is %d\n',min(tau));
fprintf('The maximum timestep is %d\n',max(tau));

%
  figure(1)
  clf ( );
  title('timesteps')
  plot ( tau, '.-' )
  print ( '-dpng', 'figure1.png' );

  figure(2)
  clf ( );
  hold ( 'on' );
    plot(ya(1),ya(2),'o')
    plot(yp(:,1),yp(:,2),'b.-')
  hold ( 'off' )
    title('phase plane')
    grid ( 'on' );
    axis ( 'equal' );
    label = sprintf ( 'Lotka-Volterra phase plot using %d intervals', iter(end) );
    title ( label );
  print ( '-dpng', 'figure2.png' );

  pause ( );

  return
end
function f_element = vector_to_element ( f_vector )

%*****************************************************************************80
%
%% vector_to_element rewrites a MATLAB vector function elementwise.
%
%  Discussion:
%
%    In MATLAB, by default, the operators '*', '/', and '^' are interpreted
%    as linear algebraic operations.  In many cases, a user simply desired
%    the elementwise operation instead, which MATLAB indicates using a 
%    dotted notation: '.*', './', and '.^'.
%
%    This function converts a string to a form that uses dotted operators.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_vector, a string presumably representing a MATLAB expression.
%
%  Output:
%
%    string f_element, a copy of f_vector in which all operators have been
%    replaced by their dotted (elementwise) form.
%
  n_vector = length ( f_vector );
  f_element = [''];

  c_prev = '';
  for iv = 1 : n_vector
    c = f_vector(iv);
    if ( ( c == '*' || c == '/' || c == '^' ) && ( c_prev ~= '.' ) )
      f_element = [ f_element, ' .', c, ' ' ];
    else
      f_element = [ f_element, c ];
    end
    c_prev = c;
  end

  return
end

