code 33302;
 procedure FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
 value N, ORDER; integer N, ORDER;
 real procedure Q, R, F;
 array X, Y, E;
 begin integer L, L1;
   real XL1, XL, H, A12, A21, B1, B2, TAU1, TAU2, CH, TL, G, YL, PP,
    E1, E2, E3, E4, E5, E6;
   array T, SUPER, SUB, CHI, GI[0:N-1];

   procedure ELEMENT MAT VEC EVALUATION 1;
   begin own real Q2, R2, F2;
    real Q1, R1, F1, H2, S12;
    if L=1 then 
    begin Q2:= Q(XL1); R2:= R(XL1); F2:= F(XL1) end;
    H2:= H/2; S12:= - 1/H;
    Q1:= Q2; Q2:= Q(XL);
    R1:= R2; R2:= R(XL);
    F1:= F2; F2:= F(XL);
    B1:= H2*F1; B2:= H2*F2;
    TAU1:= H2*R1; TAU2:= H2*R2;
    A12:= S12 + Q1/2; A21:= S12 - Q2/2
   end ELEMENT MAT VEC EV.;

   procedure ELEMENT MAT VEC EVALUATION 2;
   begin own real Q3, R3, F3;
    real Q1, Q2, R1, R2, F1, F2, S12, S13, S22, X2, H6, H15,
      C12, C32, A13, A31, A22, A23, A32, B3, TAU3;
    if L=1 then 
    begin Q3:= Q(XL1); R3:= R(XL1); F3:= F(XL1) end;

    X2:= (XL1 + XL)/2; H6:= H/6; H15:= H/1.5;
    Q1:= Q3; Q2:= Q(X2); Q3:= Q(XL);
    R1:= R3; R2:= R(X2); R3:= R(XL);
    F1:= F3; F2:= F(X2); F3:= F(XL);
    B1:= H6*F1; B2:= H15*F2; B3:= H6*F3;
    TAU1:= H6*R1; TAU2:= H15*R2; TAU3:= H6*R3;
    S12:= - 1/H/0.375; S13:= - S12/8; S22:= - 2*S12;
    A12:= S12 + Q1/1.5; A13:= S13 - Q1/6;
    A21:= S12 - Q2/1.5; A23:= S12 + Q2/1.5; A22:= S22 + TAU2;
    A31:= S13 + Q3/6; A32:= S12 - Q3/1.5;
    comment STATIC CONDENSATION;
    C12:= - A12/A22; C32:= - A32/A22;
    A12:= A13 + C12*A23; A21:= A31 + C32*A21;
    B1:= B1 + C12*B2; B2:= B3 + C32*B2;
    TAU1:= TAU1 + C12*TAU2; TAU2:= TAU3 + C32*TAU2
   end ELEMENT MAT VEC EVALUATION 2
   procedure ELEMENT MAT VEC EVALUATION 3;
   begin own real Q4, R4, F4;
    real Q1, Q2, Q3, R1, R2, R3, F1, F2, F3,
     S12, S13, S14, S22, S23, X2, X3, H12, H24,
     DET, C12, C13, C42, C43, A13, A14, A22, A23,
     A24, A31, A32, A33, A34, A41, A42, A43,
     B3, B4, TAU3, TAU4;

    if L=1 then 
    begin Q4:= Q(XL1); R4:= R(XL1); F4:= F(XL1) end;
    X2:= XL1 + 0.27639320225*H; X3:= XL - X2 + XL1;
    H12:= H/12; H24:= H/2.4;
    Q1:= Q4; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(XL);
    R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
    F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
    S12:= -4.8784183052080/H; S13:= 0.7117516385414/H;
    S14:= -.16666666666667/H; S23:= 25*S14; S22:= -2*S23;
    B1:= H12*F1; B2:= H24*F2; B3:= H24*F3; B4:= H12*F4;
    TAU1:= H12*R1; TAU2:= H24*R2; TAU3:= H24*R3; TAU4:= H12*R4;
    A12:= S12 + 0.67418082864578*Q1;
    A13:= S13 - 0.25751416197912*Q1;
    A14:= S14 + Q1/12;
    A21:= S12 - 0.67418082864578*Q2;
    A22:= S22 + TAU2;
    A23:= S23 + 0.93169499062490*Q2;
    A24:= S13 - 0.25751416197912*Q2;
    A31:= S13 + 0.25751416197912*Q3;
    A32:= S23 - 0.93169499062490*Q3;
    A33:= S22 + TAU3;
    A34:= S12 + 0.67418082864578*Q3;
    A41:= S14 - Q4/12;
    A42:= S13 + 0.25751416197912*Q4;
    A43:= S12 - 0.67418082864578*Q4;
    comment STATIC CONDENSATION;
    DET:= A22*A33 - A23*A32;
    C12:= (A13*A32 - A12*A33)/DET;
    C13:= (A12*A23 - A13*A22)/DET;
    C42:= (A32*A43 - A42*A33)/DET;
    C43:= (A42*A23 - A43*A22)/DET;
    TAU1:= TAU1 + C12*TAU2 + C13*TAU3 ;
    TAU2:= TAU4 + C42*TAU2 + C43*TAU3;
    A12:= A14 + C12*A24 + C13*A34;
    A21:= A41 + C42*A21 + C43*A31;
    B1:= B1 + C12*B2 + C13*B3;
    B2:= B4 + C42*B2 + C43*B3
   end ELEMENT MAT VEC EVALUATION 3
   procedure BOUNDARY CONDITIONS;
   if L=1 and E2 = 0 then 
   begin TAU1:= 1; B1:= E3/E1; A12:= 0 end 
   else if L=1 and E2 ^= 0 then 
   begin TAU1:= TAU1 - E1/E2; B1:= B1 - E3/E2
   end else if L=N and E5 = 0 then 
   begin TAU2:= 1; A21:= 0; B2:= E6/E4;
   end else if L=N and E5 ^= 0 then 
   begin TAU2:= TAU2 + E4/E5; B2:= B2 + E6/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;

   L:= 0; XL:= X[0];
   E1:= E[1]; E2:= E[2]; E3:= E[3]; E4:= E[4]; E5:= E[5]; E6:= E[6];
   comment ELEMENTWISE ASSEMBLAGE OF MATRIX AND VECTOR
    COMBINED WITH FORWARD BABUSKA SUBSTITUTION;
   for L:= L + 1 while L <= N do 
   begin XL1:= XL; L1:= L - 1; XL:= X[L]; H:= XL - XL1;
    if ORDER = 2 then ELEMENT MAT VEC EVALUATION 1 else 
    if ORDER = 4 then ELEMENT MAT VEC EVALUATION 2 else 
    ELEMENT MAT VEC EVALUATION 3;
    if L=1 or L=N then BOUNDARY CONDITIONS;
    FORWARD BABUSKA
   end;
   BACKWARD BABUSKA;
 end FEM LAGR;
     eop