function [x4,x5] = rkf45rt3_step(xc,xd,t,h)
%
%  This function computes a single step of Runge-Kutta-Fehlberg 4/5.
%  It returns both estimates of x(k+1).
%
%         [x4,x5] = rkf45rt3_step(xc,xd,t,h)
%
%  Input parameters:
%  -----------------
%
%     xc  :=  Value of the continuous state vector at time t.
%
%     xd  :=  Value of the discrete state vector at time t.
%
%     t   :=  Simulation time at the beginning of the step.
%
%     h   :=  Desired step length, can be either positive or negative.
%
%  Output parameters:
%  ------------------
%
%     x4  :=  4th order approximation of the state vector at time (t+h).
%
%     x5  :=  5th order approximation of the state vector at time (t+h).
%
%  Function calls:
%  ---------------
%
%  This function calls various times the function:
%
%           xcdot = cst_eq3(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.
%
xcdot   = cst_eq3(xc,xd,t);
%
xP1    = xc + (1/4)*h*xcdot;
tP1    = t + (1/4)*h;
xP1dot = cst_eq3(xP1,xd,tP1);
%
xP2    = xc + (3/32)*h*xcdot + (9/32)*h*xP1dot;
tP2    = t + (3/8)*h;
xP2dot = cst_eq3(xP2,xd,tP2);
%
xP3    = xc + (1932/2197)*h*xcdot - (7200/2197)*h*xP1dot + (7296/2197)*h*xP2dot;
tP3    = t + (12/13)*h;
xP3dot = cst_eq3(xP3,xd,tP3);
%
xP4    = xc + (439/216)*h*xcdot - 8*h*xP1dot + (3680/513)*h*xP2dot - (845/4104)*h*xP3dot;
tP4    = t + h;
xP4dot = cst_eq3(xP4,xd,tP4);
%
xP5    = xc - (8/27)*h*xcdot + 2*h*xP1dot - (3544/2565)*h*xP2dot + (1859/4104)*h*xP3dot ...
         - (11/40)*h*xP4dot;
tP5    = t + (1/2)*h;
xP5dot = cst_eq3(xP5,xd,tP5);
%
x4 = xc + (25/216)*h*xcdot + (1408/2565)*h*xP2dot + (2197/4104)*h*xP3dot - (1/5)*h*xP4dot;
x5 = xc + (16/135)*h*xcdot + (6656/12825)*h*xP2dot + (28561/56430)*h*xP3dot ...
     - (9/50)*h*xP4dot + (2/55)*h*xP5dot;
%
return