code 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 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 and E2 = 0 then 
   begin TAU1:= 1; B1:= A12:= 0 end 
   else if L=1 and E2 ^= 0 then 
   begin TAU1:= TAU1 - E1/E2
   end else if L=N and E5 = 0 then 
   begin TAU2:= 1; B2:= A21:= 0
   end else if L=N and 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 or 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;
      eop