function [xnew,xtemp,xdotnew,xdottemp] = radau_step(x,t,h)
%
%  One step of inlined Radau IIA(3) for problem [H8.2]
%  The Hessian is computed symbolically
%
R1 = 100;
R2 = 100;
R3 = 20;
C = 1e-6;
L = 0.01;
%
tol = 1e-6;
err = 1;
u3old = x(1);
iLold = x(2);
xx = [iLold;u3old;iLold;u3old];
%
while err > tol,
  iL = xx(1);
  u3 = xx(2);
  jL = xx(3);
  v3 = xx(4);
  %
  %  Inlined model equations
  %
  v0 = 10;
  u0 = 10;
  j3 = v3/R3;
  j4 = 4*v3;
  v1 = v0 - v3;
  j1 = v1/R1;
  j2 = j4 - jL;
  v2 = R2*j2;
  i3 = u3/R3;
  i4 = 4*u3;
  u1 = u0 - u3;
  i1 = u1/R1;
  i2 = i4 - iL;
  u2 = R2*i2;
  uL = u1 + u2;
  diL = uL/L;
  iC = i1 - i2 - i3;
  du3 = iC/C;
  vL = v1 + v2;
  djL = vL/L;
  jC = j1 - j2 - j3;
  dv3 = jC/C;
  iLnew = iLold + 3*h*djL/4 + h*diL/4;
  u3new = u3old + 3*h*dv3/4 + h*du3/4;
  jLnew = iLold + 5*h*djL/12 - h*diL/12;
  v3new = u3old + 5*h*dv3/12 - h*du3/12;
  u4 = u2 - u3;
  i0 = i1 + iL;
  v4 = v2 - v3;
  j0 = j1 + jL;
  %
  %  Compute the zero function for the Newton iteration
  %
  f1 = iLnew - iL;
  f2 = u3new - u3;
  f3 = jLnew - jL;
  f4 = v3new - v3;
  f = [f1;f2;f3;f4];
  %
  %  Compute symbolically the Hessian matrix for the Newton iteration
  %
  v0d1 = 0;
  u0d1 = 0;
  j3d1 = 0;
  j4d1 = 0;
  v1d1 = v0d1;
  j1d1 = v1d1/R1;
  j2d1 = j4d1;
  v2d1 = R2*j2d1;
  i3d1 = 0;
  i4d1 = 0;
  u1d1 = u0d1;
  i1d1 = u1d1/R1;
  i2d1 = i4d1 - 1;
  u2d1 = R2*i2d1;
  uLd1 = u1d1 + u2d1;
  diLd1 = uLd1/L;
  iCd1 = i1d1 - i2d1 - i3d1;
  du3d1 = iCd1/C;
  vLd1 = v1d1 + v2d1;
  djLd1 = vLd1/L;
  jCd1 = j1d1 - j2d1 - j3d1;
  dv3d1 = jCd1/C;
  iLnewd1 = 3*h*djLd1/4 + h*diLd1/4;
  u3newd1 = 3*h*dv3d1/4 + h*du3d1/4;
  jLnewd1 = 5*h*djLd1/12 - h*diLd1/12;
  v3newd1 = 5*h*dv3d1/12 - h*du3d1/12;
  %
  v0d2 = 0;
  u0d2 = 0;
  j3d2 = 0;
  j4d2 = 0;
  v1d2 = v0d2;
  j1d2 = v1d2/R1;
  j2d2 = j4d2;
  v2d2 = R2*j2d2;
  i3d2 = 1/R3;
  i4d2 = 4;
  u1d2 = u0d2 - 1;
  i1d2 = u1d2/R1;
  i2d2 = i4d2;
  u2d2 = R2*i2d2;
  uLd2 = u1d2 + u2d2;
  diLd2 = uLd2/L;
  iCd2 = i1d2 - i2d2 - i3d2;
  du3d2 = iCd2/C;
  vLd2 = v1d2 + v2d2;
  djLd2 = vLd2/L;
  jCd2 = j1d2 - j2d2 - j3d2;
  dv3d2 = jCd2/C;
  iLnewd2 = 3*h*djLd2/4 + h*diLd2/4;
  u3newd2 = 3*h*dv3d2/4 + h*du3d2/4;
  jLnewd2 = 5*h*djLd2/12 - h*diLd2/12;
  v3newd2 = 5*h*dv3d2/12 - h*du3d2/12;
  %
  v0d3 = 0;
  u0d3 = 0;
  j3d3 = 0;
  j4d3 = 0;
  v1d3 = v0d3;
  j1d3 = v1d3/R1;
  j2d3 = j4d3 - 1;
  v2d3 = R2*j2d3;
  i3d3 = 0;
  i4d3 = 0;
  u1d3 = u0d3;
  i1d3 = u1d3/R1;
  i2d3 = i4d3;
  u2d3 = R2*i2d3;
  uLd3 = u1d3 + u2d3;
  diLd3 = uLd3/L;
  iCd3 = i1d3 - i2d3 - i3d3;
  du3d3 = iCd3/C;
  vLd3 = v1d3 + v2d3;
  djLd3 = vLd3/L;
  jCd3 = j1d3 - j2d3 - j3d3;
  dv3d3 = jCd3/C;
  iLnewd3 = 3*h*djLd3/4 + h*diLd3/4;
  u3newd3 = 3*h*dv3d3/4 + h*du3d3/4;
  jLnewd3 = 5*h*djLd3/12 - h*diLd3/12;
  v3newd3 = 5*h*dv3d3/12 - h*du3d3/12;
  %
  v0d4 = 0;
  u0d4 = 0;
  j3d4 = 1/R3;
  j4d4 = 4;
  v1d4 = v0d4 - 1;
  j1d4 = v1d4/R1;
  j2d4 = j4d4;
  v2d4 = R2*j2d4;
  i3d4 = 0;
  i4d4 = 0;
  u1d4 = u0d4;
  i1d4 = u1d4/R1;
  i2d4 = i4d4;
  u2d4 = R2*i2d4;
  uLd4 = u1d4 + u2d4;
  diLd4 = uLd4/L;
  iCd4 = i1d4 - i2d4 - i3d4;
  du3d4 = iCd4/C;
  vLd4 = v1d4 + v2d4;
  djLd4 = vLd4/L;
  jCd4 = j1d4 - j2d4 - j3d4;
  dv3d4 = jCd4/C;
  iLnewd4 = 3*h*djLd4/4 + h*diLd4/4;
  u3newd4 = 3*h*dv3d4/4 + h*du3d4/4;
  jLnewd4 = 5*h*djLd4/12 - h*diLd4/12;
  v3newd4 = 5*h*dv3d4/12 - h*du3d4/12;
  %
  h11 = iLnewd1 - 1;
  h12 = iLnewd2;
  h13 = iLnewd3;
  h14 = iLnewd4;
  h21 = u3newd1;
  h22 = u3newd2 - 1;
  h23 = u3newd3;
  h24 = u3newd4;
  h31 = jLnewd1;
  h32 = jLnewd2;
  h33 = jLnewd3 - 1;
  h34 = jLnewd4;
  h41 = v3newd1;
  h42 = v3newd2;
  h43 = v3newd3;
  h44 = v3newd4 - 1;
  H = [h11,h12,h13,h14
       h21,h22,h23,h24
       h31,h32,h33,h34
       h41,h42,h43,h44];
  %
  xxnew = xx - H\f;
  err = norm(xxnew - xx);
  xx = xxnew;
end
%
xnew = [u3new;iLnew];
xtemp = [v3new;jLnew];
xdotnew = [du3;diL];
xdottemp = [dv3;djL];
%
return