function [ xp, yp ] = ode_midpoint_system ( f_string, a, b, ya, n )

%*****************************************************************************80
%
%% ode_midpoint_system is an interactive midpoint method ODE system solver.
%
%  Discussion:
%
%    This function manages an interactive computation which takes n steps
%    of the midpoint method from a to b, 
%    to estimate the solution of a system of equations y'=f(x,y).
%
%    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, 
%    and the number of steps n.
%
%    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 ( 'lotka(x,y)', 0.0, 1.0, [2.0,3.0], 10 );
%    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 number of steps: 40
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL 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. 
%
%    integer n: the number of steps to take.
%
%  Output:
%
%    real xp(n+1), yp(m,n+1): the sequence of solution estimates, which can be plotted.
%

%
%  Get f, a, b, ya, n.
%
  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 )
    n = input ( 'Enter the number of steps n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  xp = linspace ( a, b, n + 1 );
  xp = xp';
  yp = zeros ( m, n + 1 );

  dt = ( b - a ) / n;

  theta = 0.5;

  yp(:,1) = ya(:);
  for i = 1 : n
    xm = xp(i) + theta * dt; 
    ym(:,1) = yp(:,i);
    for j = 1 : 10
      ym(:,1) = yp(:,i) + theta * dt * f(xm,ym(:,1));
    end
    yp(:,i+1) = (       1.0 / theta ) * ym(:,1) ...
              + ( 1.0 - 1.0 / theta ) * yp(:,i);
  end

  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 GNU LGPL 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

