function [xnew] = bi45t_step(x,t,h,theta,tol)
%
%  This function computes a single step of Backinterpolation 4/5 with a
%  given theta value.
%
%         [xnew] = bi45t_step(x,t,h,theta,tol)
%
%  Input parameters:
%  -----------------
%
%     x     :=  Value of the state vector at time t.
%
%     t     :=  Simulation time at the beginning of the step.
%
%     h     :=  Desired step length, can be either positive or negative.
%
%     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   :=  Maximum tolerated relative error.
%
%  Output parameters:
%  ------------------
%
%     xnew  :=  Approximation of the state vector at time (t+h).
%
%  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)
%
h_left = theta*h;
h_right = (theta-1)*h;
H = hessian(x,t,h_right,4);
[x_left,dummy] = rkf45_step(x,t,h_left);
xnew = x_left;
err = 1;
while err > 0.1*tol,
  [dummy,x_right] = rkf45_step(xnew,t+h,h_right);
  xnew = xnew - H\(x_right - x_left);
  err = norm(x_right - x_left,'inf')/max([norm(x_left),norm(x_right),1.0e-10]);
end
return