%  Example of the solution of a system of ordinary differential
%  equations, including interpolation.
%  REQUIRES Rke.m, Yvalue.m, and X2rkef.m.

%  Set initial values and tolerances.
y = [1; 1];
h = 0.1;
tol = 1e-5;
threshold = [1e-5; 1e-5];

%  Define the solution output points.
nout = 11;
xout = 0:(nout-1);
    
%  Integration loop.  Rke chooses the step size for efficiency.  Yvalue
%  is used to get output by interpolation each time that Rke steps past
%  an output point xout(i), including the initial point.
x = xout(1);
nasteps = 0;                      % Initialize Rke.
yp = zeros(2,1);
for i = 1:nout
  while x <= xout(i)
    [x,y,yp,h,nasteps,flag,step,ycoeff] = ...
      Rke('X2rkef',x,y,yp,h,tol,threshold,nasteps);
    if flag ~= 0
      error(['flag = ',int2str(flag),' at x = ',num2str(x),'.']);
    end
  end

  z = Yvalue(xout(i),x,step,ycoeff);
  fprintf('Xout = %6.2f   Y(1) = %10.6f   Y(2) = %10.6f\n', ...
           xout(i),z(1),z(2));
end