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

%*****************************************************************************80
%
%% ode_trapezoidal is an interactive trapezoidal method ODE solver.
%
%  Discussion:
%
%    This function manages an interactive computation which takes n steps
%    of the trapezoidal method from a to b, to estimate the solution of y'=f(x,y).
%
%    The user enters a formula for f(x,y), the values of a and b,
%    the initial condition y(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_trapezoidal ( 'y', 0.0, 1.0, 1.0, 10 )
%    or, if called with no arguments, it will request them:
%      [ xp, yp ] = ode_trapezoidal ( );
%      Enter function formula, like x*y: y
%      Enter the starting point, a: 0.0
%      Enter the ending point, b: 1.0
%      Enter the initial condition, y(a): 1.0
%      Enter number of steps n: 10
%
%    The function is specified as a string which is either:
%    * a MATLAB expression using the arguments 'x' and 'y'
%    * the name of an M-file followed by the argument '(x,y)'.
%
%    The string should not contain any spaces between symbols, except
%    when it is passed as a function argument in quotes.
%
%    It is not necessary to use the "dot" notation for expressions involving
%    '*', '/', or '^', but it doesn't hurt either.
%
%    Examples of function specifications:
%
%      y
%      2*exp(x)-2*y
%      x*y^2
%      sin(x)*sqrt(y)/8
%      wiggle(x,y)     <-- where "wiggle.m" is a user-provided M file.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 July 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the right hand side of the differential equation.
%
%    real a, b: the starting and ending points.
%
%    real ya: the initial condition, y(a). 
%
%    integer n: the number of steps to take.
%
%  Output:
%
%    real xp(n+1), yp(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 formula, like 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 y(a): ' );
  elseif ( ischar ( ya ) )
    ya = str2num ( ya );
  end

  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 ( n + 1, 1 );

  h = ( b - a ) / n;

  yp(1) = ya;
  for i = 1 : n
%
%  Perform an iteration that we hope will solve the implicit equation
%    yp(i+1) = yp(i) + 0.5 * h * ( f(xp(i),yp(i)) + f(xp(i+1),yp(i+1)) ).
%
    z = yp(i);
    for j = 1 : 10
      z = yp(i) + 0.5 * h * ( f ( xp(i), yp(i) ) + f ( xp(i+1), z ) );
    end
    yp(i+1) = z;
  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 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

