Thefollowing are the FORTRAN programs (FORTRAN 90) thatwere used to compute the real root of a nonlinear equations system, byusing the gradients method.

…

CALL GRADIENTS (X01,X02, YT1, YT2, XP, YP, ZP, X, Y, Z, CT)

…

SUBROUTINE GRADIENTS (X01,X02, YT1, YT2, XP, YP, ZP, X, Y, Z, CT)

C   X01, X02 are the given initial values

C   YT1, YT2 are the final solutions

C   XP, YP, ZP are the coordinates of thecalculation points

C   X, Y, Z are the coordinates of the eightnodal points of the a curved quadrilateral element

C   CT is the anisotropic tensor

C   EPS is a stopping criterion

IMPLICIT NONE

DOUBLE PRECISION  X(8), Y(8), Z(8), X01, X02, YT1,YT2, XP, YP, ZP, EPS, Y1, Y2

DOUBLE PRECISION F, F1, F2, CT(3, 3)

INTEGER L

L=5000

   EPS=1.0E-35

   YT1=X01

   YT2=X02        

   CALL DSNSE(EPS, YT1,YT2, Y1, Y2, L, XP, YP, ZP, X, Y, Z, CT)

   IF(L.GT.0) THEN

   WRITE(11,*) 'iteration number = ',5000-L  

   ENDIF

   RETURN

   END

SUBROUTINE DSNSE(EPS,YT1, YT2, Y1, Y2, L, XP, YP, ZP, X, Y, Z, CT)

   IMPLICIT NONE

      DOUBLE PRECISION EPS, YT1, YT2, Y1, Y2, X(8), Y(8), Z(8), XP, YP, ZP, F, D, S, CT(3, 3)

   INTEGER L

5   CALL FS(YT1, YT2, F, Y1, Y2, XP, YP, ZP, X, Y, Z, CT)

   IF(DSQRT(F).GE.EPS) THEN

   L=L-1

   IF(L.EQ.0) RETURN

    D=Y1**2+Y2**2

     IF(D.EQ.0)THEN

    L=-1

     RETURN

     ENDIF

    S=F/D

    YT1=YT1-S*Y1

    YT2=YT2-S*Y2

     GOTO 5

     ENDIF

     RETURN

     END

SUBROUTINE FS(YT1,YT2, F, Y1, Y2, XP, YP, ZP, X, Y, Z, CT)

   IMPLICIT NONE

      DOUBLE PRECISION YT1, YT2, Y1, Y2, X(8), Y(8), Z(8), XP, YP, ZP, F, XK, YK, ZK

   DOUBLE PRECISION F1, F2, DF1, DF2, DX1, DY1, DZ1, DX2, DY2, DZ2, CT(3, 3)

   DOUBLE PRECISION Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, DF11, DF12, DF21, DF22, DX11, DY11, DZ11

   DOUBLE PRECISION DX12, DY12, DZ12, RR1, RR2, RR3

   DOUBLE PRECISION DX21, DY21, DZ21

   DOUBLE PRECISION DX22, DY22, DZ22

   DX1 =(1-YT2)*(2.D0*YT1+YT2)/4.D0*X(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*X(2)

    +           + (1+YT2)*(2.D0*YT1+YT2)/4.D0*X(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*X(4)

    +           + YT1*(YT2-1.D0)*X(5)

    +           + (1.D0-YT2**2)/2.D0*X(6)

    +           -  YT1*(1.D0+YT2)*X(7)

    +           +(YT2**2-1.D0)/2.D0*X(8)  

      DY1 = (1-YT2)*(2.D0*YT1+YT2)/4.D0*Y(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*Y(2)

    +           +(1+YT2)*(2.D0*YT1+YT2)/4.D0*Y(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*Y(4)

    +           + YT1*(YT2-1.D0)*Y(5)

    +           + (1.D0-YT2**2)/2.D0*Y(6)

    +           -  YT1*(1.D0+YT2)*Y(7)

    +           +(YT2**2-1.D0)/2.D0*Y(8)  

       DZ1 = (1-YT2)*(2.D0*YT1+YT2)/4.D0*Z(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*Z(2)

    +           +(1+YT2)*(2.D0*YT1+YT2)/4.D0*Z(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*Z(4)

    +           + YT1*(YT2-1.D0)*Z(5)

    +           + (1.D0-YT2**2)/2.D0*Z(6)

    +           -  YT1*(1.D0+YT2)*Z(7)

    +           + (YT2**2-1.D0)/2.D0*Z(8)

       DX2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*X(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*X(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*X(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*X(4)

    +           + (YT1**2-1.D0)/2.D0*X(5)

    +           - YT2*(1.D0+YT1)*X(6)

    +           + (1.D0-YT1**2)/2.D0*X(7)

    +           +YT2*(YT1-1.D0)*X(8)  

       DY2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*Y(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*Y(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*Y(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*Y(4)

    +           + (YT1**2-1.D0)/2.D0*Y(5)

    +           - YT2*(1.D0+YT1)*Y(6)

    +           + (1.D0-YT1**2)/2.D0*Y(7)

    +           + YT2*(YT1-1.D0)*Y(8)  

       DZ2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*Z(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*Z(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*Z(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*Z(4)

    +           + (YT1**2-1.D0)/2.D0*Z(5)

    +           - YT2*(1.D0+YT1)*Z(6)

    +           + (1.D0-YT1**2)/2.D0*Z(7)

    +           + YT2*(YT1-1.D0)*Z(8)

       Q1 = (1-YT1)*(1-YT2)*(-YT1-YT2-1.D0)/4.D0

       Q2 = (1+YT1)*(1-YT2)*(YT1-YT2-1.D0)/4.D0

       Q3 = (1+YT1)*(1+YT2)*(YT1+YT2-1.D0)/4.D0

       Q4 = (1-YT1)*(1+YT2)*(-YT1+YT2-1.D0)/4.D0

       Q5 = (1-YT1**2)*(1-YT2)/2.D0

       Q6 = (1+YT1)*(1-YT2**2)/2.D0

       Q7 = (1+YT2)*(1-YT1**2)/2.D0

       Q8 = (1-YT1)*(1-YT2**2)/2.D0

       XK=Q1*X(1)+Q2*X(2)+Q3*X(3)+Q4*X(4)+Q5*X(5)+Q6*X(6)+Q7*X(7)+Q8*X(8)

       YK=Q1*Y(1)+Q2*Y(2)+Q3*Y(3)+Q4*Y(4)+Q5*Y(5)+Q6*Y(6)+Q7*Y(7)+Q8*Y(8)

       ZK=Q1*Z(1)+Q2*Z(2)+Q3*Z(3)+Q4*Z(4)+Q5*Z(5)+Q6*Z(6)+Q7*Z(7)+Q8*Z(8)

        RR1=XK-XP

        RR2=YK-YP

        RR3=ZK-ZP

        F1=CT(1,1)*DX1*RR1+CT(1,2)*DX1*RR2+CT(1,3)*DX1*RR3

    +      +CT(2,1)*DY1*RR1+CT(2,2)*DY1*RR2+CT(2,3)*DY1*RR3

    +      +CT(3,1)*DZ1*RR1+CT(3,2)*DZ1*RR2+CT(3,3)*DZ1*RR3

       F2=CT(1,1)*DX2*RR1+CT(1,2)*DX2*RR2+CT(1,3)*DX2*RR3

    +     +CT(2,1)*DY2*RR1+CT(2,2)*DY2*RR2+CT(2,3)*DY2*RR3

    +     +CT(3,1)*DZ2*RR1+CT(3,2)*DZ2*RR2+CT(3,3)*DZ2*RR3

       F=F1*F1+F2*F2

       DX11=(1-YT2)/2.D0*X(1)+(1-YT2)/2.D0*X(2)+(1+YT2)/2.D0*X(3)+(1+YT2)/2.D0*X(4)+(YT2-1)*X(5)+(-YT2-1)*X(7)

       DY11=(1-YT2)/2.D0*Y(1)+(1-YT2)/2.D0*Y(2)+(1+YT2)/2.D0*Y(3)+(1+YT2)/2.D0*Y(4)+(YT2-1)*Y(5)+(-YT2-1)*Y(7)

       DZ11=(1-YT2)/2.D0*Z(1)+(1-YT2)/2.D0*Z(2)+(1+YT2)/2.D0*Z(3)+(1+YT2)/2.D0*Z(4)+(YT2-1)*Z(5)+(-YT2-1)*Z(7)

       DF11=CT(1,1)*DX11*RR1+CT(1,1)*DX1*DX1+CT(1,2)*DX11*RR2+CT(1,2)*DX1*DY1+CT(1,3)*DX11*RR3+CT(1,3)*DX1*DZ1

    +           +CT(2,1)*DY11*RR1+CT(2,1)*DY1*DX1+CT(2,2)*DY11*RR2+CT(2,2)*DY1*DY1+CT(2,3)*DY11*RR3+CT(2,3)*DY1*DZ1

    +          +CT(3,1)*DZ11*RR1+CT(3,1)*DZ1*DX1+CT(3,2)*DZ11*RR2+CT(3,2)*DZ1*DY1+CT(3,3)*DZ11*RR3+CT(3,3)*DZ1*DZ1        

       DX12=(1-2.D0*YT1-2.D0*YT2)/4.D0*X(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*X(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*X(3)

    +          +(-1+2.D0*YT1-2.D0*YT2)/4.D0*X(4)

    +           +YT1*X(5)-YT2*X(6)-YT1*X(7)+YT2*X(8)

       DY12=(1-2.D0*YT1-2.D0*YT2)/4.D0*Y(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*Y(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*Y(3)

    +           +(-1+2.D0*YT1-2.D0*YT2)/4.D0*Y(4)

    +           +YT1*Y(5)-YT2*Y(6)-YT1*Y(7)+YT2*Y(8)

       DZ12=(1-2.D0*YT1-2.D0*YT2)/4.D0*Z(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*Z(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*Z(3)

    +          +(-1+2.D0*YT1-2.D0*YT2)/4.D0*Z(4)

    +           +YT1*Z(5)-YT2*Z(6)-YT1*Z(7)+YT2*Z(8)

       DF12=CT(1,1)*DX12*RR1+CT(1,1)*DX1*DX2+CT(1,2)*DX12*RR2+CT(1,2)*DX1*DY2+CT(1,3)*DX12*RR3+CT(1,3)*DX1*DZ2    !!!

    +          +CT(2,1)*DY12*RR1+CT(2,1)*DY1*DX2+CT(2,2)*DY12*RR2+CT(2,2)*DY1*DY2+CT(2,3)*DY12*RR3+CT(2,3)*DY1*DZ2

    +          +CT(3,1)*DZ12*RR1+CT(3,1)*DZ1*DX2+CT(3,2)*DZ12*RR2+CT(3,2)*DZ1*DY2+CT(3,3)*DZ12*RR3+CT(3,3)*DZ1*DZ2    

       DX21=DX12

       DY21=DY12

       DZ21=DZ12        

       DF21=CT(1,1)*DX21*RR1+CT(1,1)*DX2*DX1+CT(1,2)*DX21*RR2+CT(1,2)*DX2*DY1+CT(1,3)*DX21*RR3+CT(1,3)*DX2*DZ1    !!!

    +           +CT(2,1)*DY21*RR1+CT(2,1)*DY2*DX1+CT(2,2)*DY21*RR2+CT(2,2)*DY2*DY1+CT(2,3)*DY21*RR3+CT(2,3)*DY2*DZ1

    +          +CT(3,1)*DZ21*RR1+CT(3,1)*DZ2*DX1+CT(3,2)*DZ21*RR2+CT(3,2)*DZ2*DY1+CT(3,3)*DZ21*RR3+CT(3,3)*DZ2*DZ1        

       DX22=(1-YT1)/2.D0*X(1)+(1+YT1)/2.D0*X(2)+(1+YT1)/2.D0*X(3)+(1-YT1)/2.D0*X(4)+(-YT1-1)*X(6)+(YT1-1)*X(8)

       DY22=(1-YT1)/2.D0*Y(1)+(1+YT1)/2.D0*Y(2)+(1+YT1)/2.D0*Y(3)+(1-YT1)/2.D0*Y(4)+(-YT1-1)*Y(6)+(YT1-1)*Y(8)

       DZ22=(1-YT1)/2.D0*Z(1)+(1+YT1)/2.D0*Z(2)+(1+YT1)/2.D0*Z(3)+(1-YT1)/2.D0*Z(4)+(-YT1-1)*Z(6)+(YT1-1)*Z(8)

       DF22=CT(1,1)*DX22*RR1+CT(1,1)*DX2*DX2+CT(1,2)*DX22*RR2+CT(1,2)*DX2*DY2+CT(1,3)*DX22*RR3+CT(1,3)*DX2*DZ2    !!!

    +          +CT(2,1)*DY22*RR1+CT(2,1)*DY2*DX2+CT(2,2)*DY22*RR2+CT(2,2)*DY2*DY2+CT(2,3)*DY22*RR3+CT(2,3)*DY2*DZ2

    +          +CT(3,1)*DZ22*RR1+CT(3,1)*DZ2*DX2+CT(3,2)*DZ22*RR2+CT(3,2)*DZ2*DY2+CT(3,3)*DZ22*RR3+CT(3,3)*DZ2*DZ2

       Y1=2.D0*(F1*DF11+F2*DF21)

       Y2=2.D0*(F1*DF12+F2*DF22)

     RETURN

     END