function [tout, yout] = ode23(FunFcn, t0, tfinal, y0, tol, trace)
%
% kris' version for the cannon problem - need to quit when second
%       component goes negative
%
%ODE23	Integrate a system of ordinary differential equations using
%	2nd and 3rd order Runge-Kutta formulas.  See also ODE45 and
%	ODEDEMO.M.
%	[T,Y] = ODE23('yprime', T0, Tfinal, Y0) integrates the system
%	of ordinary differential equations described by the M-file
%	YPRIME.M over the interval T0 to Tf and using initial
%	conditions Y0.
%	[T, Y] = ODE23(F, T0, Tfinal, Y0, TOL, 1) uses tolerance TOL
%	and displays status while the integration proceeds.
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
% t0    - Initial value of t.
% tfinal- Final value of t.
% y0    - Initial value column-vector.
% tol   - The desired accuracy. (Default: tol = 1.e-3).
% trace - If nonzero, each step is printed. (Default: trace = 0).
%
% OUTPUT:
% T  - Returned integration time points (column-vector).
% Y  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).

%   C.B. Moler, 3-25-87.
%   Copyright (c) 1987 by the MathWorks, Inc.
%   All rights reserved.

% Initialization
pow = 1/3;
if nargin < 5, tol = 0.001; end
if nargin < 6, trace = 0; end

% Initialization
t = t0;
hmax = (tfinal - t)/5;
hmin = (tfinal - t)/20000;
h = (tfinal - t)/100;
y = y0(:);
tout = t;
yout = y.';
tau = tol * max(norm(y, 'inf'), 1);

if trace
   clc, t, h, y
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      s1 = feval(FunFcn, t, y);
      s2 = feval(FunFcn, t+h, y+h*s1);
      s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4);

      % Estimate the error and the acceptable error
      delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
      tau = tol*max(norm(y,'inf'),1.0);

     % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         y = y + h*(s1 + 4*s3 + s2)/6;
         tout = [tout; t];
         yout = [yout; y.'];
      end
      if trace
         home, t, h, y
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.9*h*(tau/delta)^pow);
      end

% kscheck if second component is negative
   if y(2) < 0 , break, end;

   end; %   while (t < tfinal) & (h >= hmin)

%ks   if (t < tfinal)
%ks      disp('SINGULARITY LIKELY.')
%ks      t
%ks   end