%  Example of the solution of a system of ordinary differential equations.
%  REQUIRES Rke.m and X1rkef.m.

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

%  Integration loop.  Rke chooses the step size for efficiency until it is
%  about to step past xout, at which time the step size is reduced to produce
%  an answer at xout.  Because the answer is not computed by interpolation, the
%  optional output arguments step and ycoeff can be omitted from the call list.
nasteps = 0;
yp = zeros(2,1);
while x < xout
  if (x+h) > xout
    h = xout - x;
  end
  [x,y,yp,h,nasteps,flag] = Rke('X1rkef',x,y,yp,h,tol,threshold,nasteps);
  if flag ~= 0
    error(['flag = ',int2str(flag),' at x = ',num2str(x),'.']);
  end
end

%  Display the results.
fprintf('At xout = %f\n',xout);
fprintf('The numerical solution is %f %f\n',y(1),y(2));
fprintf('The true solution is      %f %f\n',exp(xout),exp(-xout));