procedure FIND GRADIENTS AT NODES(X,Y,GRAD,LBX,UBX);
   value LBX,UBX;   integer LBX,UBX;
     array X,Y,GRAD;
     comment X AND Y ARE ARRAYS WITH SUBSCRIPTS FROM LBX
     TO UBX GIVING THE ABSCISSAE AND ORDINATES RESPECT-
     IVELY OF THE DATA POINTS.  THESE SHOULD BE IN
     ASCENDING ORDER OF ABSCISSAE, WITH NO TWO ABSCISSAE
     EQUAL.  THE ARRAYS MUST CONTAIN AT LEAST FOUR
     POINTS.
     GRAD IS AN ARRAY HAVING THE SAME SUBSCRIPTS AS X
     AND Y INTO WHICH WILL BE PLACED THE CALCULATED
     GRADIENTS OF THE REQUIRED CUBICS AT EACH OF THE
     DATA POINTS.
     LBX AND UBX ARE INTEGERS GIVING THE MINIMUM AND
     MAXIMUM VALUES RESPECTIVELY OF THE SUBSCRIPTS OF
     THE ARRAYS X, Y, AND GRAD;
begin
   integer I,ILESS2,ILESS1,IPLUS1,IPLUS2;
   real X0,X1,X2,X3,X4,Y2,PROD1,PROD2,NUM,DENOM,G,
     COEFF2,XDIFF,XPROD,WEIGHT;
   for I ≔ LBX step 1 until UBX do
      begin
         comment SPECIAL TREATMENT IS NEEDED AT END POINT;
         ILESS1 ≔ if I>LBX then I-1 else I+3;
         IPLUS1 ≔ if I<UBX then I+1 else I-3;
         X2 ≔ X[I];    Y2 ≔ Y[I];
         X1 ≔ X[ILESS1]-X2;    X3 ≔ X[IPLUS1]-X2;
         comment FIRST FIT A QUADRATIC THROUGH X1,X2,X3;
         PROD1 ≔ (Y[ILESS1]-Y2)×X3;
         PROD2 ≔ (Y[IPLUS1]-Y2)×X1;
         DENOM ≔ X1×X3×(X[ILESS1]-X[IPLUS1]);
         G ≔ (X1×PROD2-X3×PROD1)/DENOM;
         COEFF2 ≔ (PROD1-PROD2)/DENOM;
         comment IF X0 EXISTS, FIND ITS CONTRIBUTION TO THE
           CUBIC ADJUSTMENT;
         if I⩽LBX+1 then NUM ≔ DENOM ≔ 0·0 else
         begin
            ILESS2 ≔ I-2;    X0 ≔ X[ILESS2]-X2;
            XDIFF ≔ X[ILESS2]-X[ILESS1];
            XPROD ≔ X0×XDIFF×(X[ILESS2]-X[IPLUS1]);
            WEIGHT ≔ XPROD/(XDIFF×XDIFF);
            NUM ≔ WEIGHT×(Y[ILESS2]-Y2-X0×(G+X0×COEFF2));
            DENOM ≔ WEIGHT×XPROD
         end;
         comment IF X4 EXISTS, FIND ITS CONTRIBUTION TO THE
           CUBIC ADJUSTMENT;
         if I<UBX-1 then
         begin
            IPLUS2 ≔ I+2;    X4 ≔ X[IPLUS2]-X2;
            XDIFF ≔ X[IPLUS2]-X[IPLUS1];
            XPROD ≔ X4×XDIFF×(X[IPLUS2]-X[ILESS1]);
            WEIGHT ≔ XPROD/(XDIFF×XDIFF);
            NUM ≔ NUM+WEIGHT×(Y[IPLUS2]-Y2-X4×(G+X4×COEFF2));
            DENOM ≔ DENOM+WEIGHT×XPROD
         end;
         GRAD[I] ≔ G+NUM×X1×X3/DENOM
      end
end OF PROCEDURE FIND GRADIENTS AT NODES;
procedure COEFFS(X,Y,GRAD,I,C0,C1,C2,C3);
   value I;   integer I;
     real C0,C1,C2,C3;
     array X,Y,GRAD;
     comment THIS PROCEDURE CALCULATES THE COEFFICIENTS
     OF THE CUBIC (IN THE RANGE X[I] TO X[I+1]) WHICH
     HAS VALUES Y[I], Y[I 1] AND GRADIENTS GRAD[I],
     GRAD[I+1] AT X[I], X[I+1] RESPECTIVELY.
     THE CUBIC TAKES THE FORM
     C0 + C1*(X-X[I]) + C2*(X-X[I])'2 + C3*(X-X[I])'3;
begin
   real H,DY;
   C0 ≔ Y[I];
   H ≔ X[I+1]-X[I];    DY ≔ Y[I+1]-C0;
   C1 ≔ GRAD[I];
   C2 ≔ (3·0×DY-H×(2·0×C1+GRAD[I+1]))/(H×H);
   C3 ≔ (H×(C1+GRAD[I+1])-2·0×DY)/H;
end OF PROCEDURE COEFFS;