function [y,xc,xd,tout] = rkf45rt(xc0,xd0,t,tol)
%
%  This function simulates a discontinuous system using
%  Runge-Kutta-Fehlberg 4/5 with a root solver and
%  Gustaffsson step-size control.
%
%         [y,xc,xd,tout] = rkf45rt(xc0,xd0,t,tol)
%
%  Input parameters:
%  -----------------
%
%     xc0  :=  Initial value of the continuous state vector at time t.
%
%     xd0  :=  Initial value of the discrete state vector at time t.
%
%     t    :=  Row vector of communication time values.
%
%     tol  :=  Desired tolerance.
%
%  Output parameters:
%  ------------------
%
%     y    :=  Matrix of output values at communication instances.
%              Each row denotes one output variable.
%              Each column denotes one communication instance.
%
%     xc   :=  Matrix of continuous state values at communication
%              instances.
%
%     xd   :=  Matrix of discrete state values at communication
%              instances.
%
%     tout :=  Row vector of communication instances.
%              These are the same time instances as those given in t
%              augmented by the event time instances.
%              Each event time instance is contained twice in the tout
%              vector, as the output and state variable values are
%              logged immediately before and immediately after each
%              event.
%
%  User-coded functions:
%  ---------------------
%
%           xcdot = cst_eq(xc,xd,t)
%
%  to compute the state derivative at a given time, t, for a given
%  value of the continuous state vector, xc, and for a given value of
%  the discrete state vector, xd, using the state-space model of the
%  system to be simulated.
%
%           xdnew = dst_eq(xc,xd,t,evt_nbr)
%
%  to compute the new discrete state, xdnew, at a given event time, t,
%  for a given value of the continuous state vector, xc, and for a
%  given value of the discrete state vector, xd.  The last input,
%  evt_nbr, determines the type of event.  If evt_nbr > 0, we are
%  dealing with a time event, and evt_nbr is the event number.
%  If evt_nbr < 0, we are dealing with a state event, and -evt_nbr
%  is the index of the zero-crossing function that triggered the event.
%
%            f = zcf(xc,xd,t)
%
%  to compute the values of the zero-crossing functions, f, at a given
%  time, t, and for a given state, [xc,xd].
%
%             y = out_eq3(xc,xd,t)
%
%  to compute the values of the output variables, y, at a given
%  time, t, and for a given state, [xc,xd].
%
%  Internal functions:
%  -------------------
%
%             [xc4,xc5] = rkf45rt_step(xc,xd,t,h)
%
%  computes one step of Runge-Kutta-Fehlberg 4/5 starting at time t and
%  continuous state xc with step size h.  The discrete state, xd, is
%  treated like a parameter vector.
%
%             tnew = reg_falsi([t1,t2],[f1,f2])
%
%  computes one iteration step of regula falsi given the left and right
%  time values, t1 and t2, and the left and right zero-crossing
%  function vectors, f1 and f2.  tnew is the next evaluation time.
%
%             iter = zc_iter([f1,f2],tol)
%
%  checks whether one or several zero crossings have taken place.
%  The routine returns iter, an index with the following meaning:
%
%    iter = 0 :  no zero crossing took place, no iteration needed
%         = 1 :  one or more zero crossings took place, we need
%                iteration
%         < 0 :  one zero crossing took place, and we have
%                convergence.  -iter is the index of the zero crossing
%                function that triggered the state event
%
%             [] = push_evt(te,evt_nbr)
%
%  places a time event of type evt_nbr to occur at time te on the
%  event calendar
%
%             [te,evt_nbr] = pull_evt()
%
%  removes the first time event from the event calendar and returns
%  its time of occurrence, te, and the event number, evt_nbr
%
%             [te,evt_nbr] = query_evt()
%
%  returns the information concerning the first time event without
%  removing the event from the event calendar
%
%  Set the initial state
%
t1 = 0;
xc1 = xc0;
xd1 = xd0;
y1 = out_eq3(xc1,xd1,t1);
f1 = zcf(xc1,xd1,t1);
%
%  Log the initial state
%
tout = t1;
xc = xc1;
xd = xd1;
y = y1;
%
%  Determine:
%    the first communication time, tc,
%    the first time event occurrence, te,
%    the final time, tf, and
%    the first scheduled output time, to
%
i = 1;
tc = t(i);
if tc == t1,
  i = 2;
  tc = t(i);
end
n = length(t);
tf = t(n);
dt = tf - t1;
tolt = tol*dt;
te = query_evt();
if te > 0,
  to = min([tc,te]);
 else
  to = tc;
end
%
%  Determine the first step size, h
%
h = min([to - t1,dt/100]);
if t1 + 1.05*h > to,
  h = to - t1;
end
errl = 0;
%
%  Simulate using RKF4/5
%
while t1 < tf,
  [xc4,xc5] = rkf45rt_step(xc1,xd1,t1,h);
  err = norm(xc4 - xc5,'inf')/max([norm(xc4),norm(xc5),1.0e-10]);
  %
  %  Check first, if step needs to be rejected due to lack of accuracy
  %
  if err > tol,
    h = (0.8*tol/err)^(0.2) * h;
    errl = 0;
   else
    %
    %  Accept the current step temporarily
    %
    t2 = t1 + h;
    xc2 = xc5;
    %
    %  Check next whether we need to iterate on a zero-crossing 
    %
    f2 = zcf(xc2,xd1,t2);
    iter = zc_iter([f1,f2],tol);
    if iter == 1,
      %
      %  We need to reject the step and iterate
      %
      tnew = reg_falsi([t1,t2],[f1,f2]);
      h = tnew - t1;
      errl = 0;
     else
      %
      %  Accept the current step permanently
      %
      t1 = t2;
      xc1 = xc2;
      f1 = f2;
      %
      %  Check whether a time event needs to be handled
      %
      if abs(te - t1) < tolt,
        y1 = out_eq3(xc1,xd1,t1);
        tout = [tout,t1];
        y = [y,y1];
        xc = [xc,xc1];
        xd = [xd,xd1];
        [te,evt_nbr] = pull_evt();
        xd1 = dst_eq(xc1,xd1,t1,evt_nbr);
        te = query_evt();
        tout = [tout,t1];
        y = [y,y1];
        xc = [xc,xc1];
        xd = [xd,xd1];
      end,
      %
      %  Check whether a state event needs to be handled
      %
      if iter < 0,
        y1 = out_eq3(xc1,xd1,t1);
        tout = [tout,t1];
        y = [y,y1];
        xc = [xc,xc1];
        xd = [xd,xd1];
        xd1 = dst_eq(xc1,xd1,t1,iter);
        iter = 0;
        tout = [tout,t1];
        y = [y,y1];
        xc = [xc,xc1];
        xd = [xd,xd1];
      end,
      %
      %  Check whether it is time to log the output
      %
      if abs(tc - t1) < tolt,
        y1 = out_eq3(xc1,xd1,t1);
        tout = [tout,t1];
        y = [y,y1];
        xc = [xc,xc1];
        xd = [xd,xd1];
        if i < n,
          i = i + 1;
          tc = t(i);
        end,
      end,
      %
      %  Compute the new step size using Gustaffsson step-size
      %  control
      %
      if errl > 0,
        h = (0.8*tol/err)^(0.06) * (errl/err)^(0.08) * h;
      else
        h = (0.8*tol/err)^(0.2) * h;
      end,
      errl = err;
      if te > 0,
        to = min([tc,te]);
       else
        to = tc;
      end,
      if t1 + 1.05*h > to,
        h = to - t1;
        errl = 0;
      end,
    end,
  end,
end
%
return