comment ================== 33314 ================= ;
procedure NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
   integer N, NC;
     real procedure F, FY, FZ;
     array X, Y, E;
begin integer L, L1, IT;
   real XL1, XL, H, A12, A21, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
     PLM, PRM, PL1, PL3, PL1PL2, PL1PL3, PL2PL2, PL2PL3,
     PR1PR2, PR1PR3, PR2PR3, PL1QL2, PL1QL3, PL2QL1, PL2QL2, PL2QL3,
     PL3QL1, PL3QL2, PR1QR2, PR1QR3, PR2QR1, PR2QR2, PR2QR3, PR3QR1,
     PR3QR2, H2RM, ZL1, ZL, E1, E2, E3, E4, E5, E6, EPS, RHO;
   array T, SUPER, SUB, CHI, GI[0:N-1], Z[0:N];

   procedure DUPVEC(L, U, S, A, B); code 31030;

   procedure ELEMENT MAT VEC EVALUATION 1;
   begin real  XM, VL, VR, WL, WR, PR, QM, RM, FM, XL12, XL1XL, XL2, ZM, ZACCM;
      if NC = 0 then VL ≔ VR ≔ 0·5 else if NC = 1 then
      begin VL ≔ (XL1 × 2 + XL)/6; VR ≔ (XL1 + XL × 2)/6 end else
      begin XL12 ≔ XL1 × XL1/12; XL1XL ≔ XL1 × XL/6; XL2 ≔ XL × XL/12;
         VL ≔ 3 × XL12 + XL1XL + XL2;
         VR ≔ 3 × XL2 + XL1XL + XL12
      end;
      WL ≔ H × VL; WR ≔ H × VR; PR ≔ VR/(VL + VR);
      XM ≔ XL1 + H × PR; ZM ≔ PR × ZL + (1 - PR) × ZL1;
      ZACCM ≔ (ZL - ZL1)/H ; QM ≔ FZ(XM, ZM, ZACCM);
      RM ≔ FY(XM, ZM, ZACCM); FM ≔ F(XM, ZM, ZACCM);
      TAU1 ≔ WL × RM; TAU2 ≔ WR × RM;
      B1 ≔ WL × FM - ZACCM × (VL + VR); B2 ≔ WR × FM + ZACCM × (VL + VR);
      A12 ≔ - (VL + VR)/H + VL × QM + (1 - PR) × PR × RM × (WL + WR);
      A21 ≔ - (VL + VR)/H - VR × QM + (1 - PR) × PR × RM × (WL + WR);
   end ELEM. M.V. EV.;
   procedure BOUNDARY CONDITIONS;
      if L = 1 ∧ E2 = 0 then
      begin TAU1 ≔ 1; B1 ≔ A12 ≔ 0 end
        else if L = 1 ∧ E2 ≠ 0 then
      begin TAU1 ≔ TAU1 - E1/E2
      end else if L = N ∧ E5 = 0 then
      begin TAU2 ≔ 1; B2 ≔ A21 ≔ 0
      end else if L = N ∧ E5 ≠ 0 then
      begin TAU2 ≔ TAU2 + E4/E5
      end B.C.1;

   procedure FORWARD BABUSKA;
      if L = 1 then
      begin CHI[0] ≔ CH ≔ TL ≔ TAU1; T[0] ≔ TL;
         GI[0] ≔ G ≔ YL ≔ B1; Y[0] ≔ YL;
         SUB[0] ≔ A21; SUPER[0] ≔ A12;
         PP ≔ A21/(CH - A12); CH ≔ TAU2 - CH × PP;
         G ≔ B2 - G × PP; TL ≔ TAU2; YL ≔ B2
      end else
      begin CHI[L1] ≔ CH ≔ CH + TAU1;
         GI[L1] ≔ G ≔ G + B1;
         SUB[L1] ≔ A21; SUPER[L1] ≔ A12;
         PP ≔ A21/(CH - A12); CH ≔ TAU2 - CH × PP;
         G ≔ B2 - G × PP; T[L1] ≔ TL + TAU1; TL ≔ TAU2;
         Y[L1] ≔ YL +  B1; YL ≔ B2
      end FORWARD BABUSKA;

   procedure BACKWARD BABUSKA;
   begin PP ≔ YL; Y[N] ≔ G/CH;
      G ≔ PP; CH ≔ TL; L ≔ N;
      for L ≔ L - 1 while L ⩾ 0 do
         begin PP ≔ SUPER[L]/(CH - SUB[L]);
            TL ≔ T[L]; CH ≔ TL - CH × PP;
            YL ≔ Y[L]; G ≔ YL - G × PP;
            Y[L] ≔ (GI[L] + G - YL)/(CHI[L] + CH - TL)  ;
         end
   end BACKWARD BABUSKA;

   DUPVEC(0, N, 0, Z, Y);
   E1 ≔ E[1]; E2 ≔ E[2]; E3 ≔ E[3]; E4 ≔ E[4]; E5 ≔ E[5]; E6 ≔ E[6];
   for IT ≔ 1, IT + 1 while EPS > RHO do
      begin L ≔ 0; XL ≔ X[0]; ZL ≔ Z[0];
         for L ≔ L + 1 while L ⩽ N do
            begin XL1 ≔ XL; L1 ≔ L - 1; XL ≔ X[L]; H ≔ XL - XL1;
               ZL1 ≔ ZL; ZL ≔ Z[L];
               ELEMENT MAT VEC EVALUATION 1;
               if L = 1 ∨ L = N then BOUNDARY CONDITIONS;
               FORWARD BABUSKA
            end;
         BACKWARD BABUSKA;
         EPS ≔ 0; RHO ≔ 1;
         for L ≔ 0 step 1 until N do
            begin RHO ≔ RHO + ABS(Z[L]);
               EPS ≔ EPS + ABS(Y[L]); Z[L] ≔ Z[L] - Y[L]
            end;
         RHO ≔ ⏨-14 × RHO
      end;
   DUPVEC(0, N, 0, Y, Z)
end NONLIN FEM LAG SKEW;