function [xnew,err,tnew,hnew] = bi45tv_step(x,t,h,errl,theta,tol)
%
%  This function computes a single step of Backinterpolation 4/5 with a
%  given theta value and with step size control.
%
%         [xnew,err,tnew,hnew,nfct] = bi45tv_step(x,t,h,errl,theta,tol)
%
%  Input parameters:
%  -----------------
%
%     x     :=  Value of the state vector at time t.
%
%     t     :=  Simulation time at the beginning of the step.
%
%     h     :=  Current step length, can be either positive or negative.
%
%     errl  :=  Relative error of the previous step.  Set equal zero
%               during the first step.
%
%     theta :=  Theta value to be used in the method.  The forward semistep
%               integrates across theta*h time units, whereas the subsequent
%               backward semistep integrates across the remaining (1-theta)*h
%               time units.
%
%     tol   :=  Largest tolerated relative error.
%
%  Output parameters:
%  ------------------
%
%     xnew  :=  Approximation of the state vector at time (t+h).
%
%     err   :=  Relative error of the current step.  If the step was
%               rejected, err is returned with a negative sign.
%
%     tnew  :=  Time at the end of the step, equal to t+h except if
%               the step was repeated with a smaller step size.
%
%     hnew  :=  Step size proposed for use during the next step.
%
%  Function calls:
%  ---------------
%
%  This function calls the function:
%
%           [x4,x5] = rkf45_step(xx,tt,hh)
%
%  to compute the forward and backward semisteps composing the method.
%  It uses the 4th order approximation, x4, during the forward semistep,
%  and the 5th order approximation, x5, during the backward semistep.
%
%  It also calls the function:
%
%           [H] = hessian(xx,tt,hh,order)
%
%  to compute the Hessian of the system (order=5).
%
%  Indirect calls are to the functions:
%
%           [xxdot] = st_eq(xx,tt)
%           [J] = jacobian(xx,tt)
%
err = 2*tol;
rep = 0;
while err > tol,
  h_left = theta*h;
  h_right = (theta-1)*h;
  H = hessian(x,t,h_right,5);
  [x_left4,x_left5] = rkf45_step(x,t,h_left);
  err = norm(x_left4 - x_left5,'inf')/max([norm(x_left4),norm(x_left5),1.0e-10]);
  if err > tol,
    h = (0.8*tol/err)^(0.2) * h;
    rep = 1;
   else
    if errl > 0  &  rep == 0,
      hnew = (0.8*tol/err)^(0.06) * (errl/err)^(0.08) * h;
     else
      hnew = (0.8*tol/err)^(0.2) * h;
    end,
  end,
end
if rep == 1,
  err = -err;
end
tnew = t + h;
xnew = x_left4;
err2 = 1;
while err2 > 0.1*tol,
  [dummy,x_right] = rkf45_step(xnew,t+h,h_right);
  xnew = xnew - H\(x_right - x_left4);
  err2 = norm(x_right - x_left4,'inf')/max([norm(x_left4),norm(x_right),1.0e-10]);
end
return