function [xnew,hnew] = radau_stepv(x,t,h)
%
%  One step of inlined Radau IIA(3) for problem [H8.3]
%  with Fehlberg (optimistic) step-size control
%  The Hessian is computed symbolically
%
%  One step of of length h step-size controlled Radau IIA(3)
%  consists of two semi-steps of length h/2 of Radau IIA(3)
%  without step-size control
%
h2 = h/2;
tol = 1e-6;
delta = 1e-10;
%
%  Compute the two semi-steps:
%
x1 = x;
[x3,x2] = radau_step(x,t,h2);
[xnew,x4,xnewdot,x4dot] = radau_step(x3,t+h2,h2);
%
%  Compute the error method
%
xblend = (-x1 + 2*x2 + 14*x3 - 2*x4 + (11*x4dot + 3*xnewdot)*h2)/13;
%
%  Compute the new step size
%
hnew = (tol*max([norm(xnew),norm(xblend),delta])/norm(xnew - xblend))^(1/3)*h;
%
return