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

%*****************************************************************************80
%
%% ode_euler_system() is an interactive Euler-method ODE system solver.
%
%  Discussion:
%
%    This function manages an interactive computation which takes n steps
%    of Euler's method from a to b, to estimate the solution of 
%    a system of differential 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_euler_system ( 'lotka(x,y)', 0.0, 50.0, [4,2], 8000 )
%    or, if called with no arguments, it will request them:
%      [ xp, yp ] = ode_euler_system ( );
%      Enter function formula, like lotka(x,y): lotka(x,y)
%      Enter the starting point, a: 0.0
%      Enter the ending point, b: 50.0
%      Enter the initial condition, y(a), such as [1,2,3]: [4,2]
%      Enter number of steps n: 8000
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 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 condition vector, y(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 = str2func ( [ '@(x,y)', f_string ] );

  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 );

  h = ( b - a ) / n;

  yp(:,1) = ya(:);
  for i = 1 : n
    yp(:,i+1) = yp(:,i) + h * f ( xp(i), yp(:,i) );
  end

  return
end

