SUBROUTINE DERPAR(N, X, XLOW, XUPP, EPS, W, INITAL, ITIN, HH,     DER   10
     * HMAX, PREF, NDIR, E, MXADMS, NCORR, NCRAD, NOUT, OUT,
     * MAXOUT, NPRNT)
      DIMENSION X(11), XLOW(11), XUPP(11), W(11), HMAX(11),
     * PREF(11), NDIR(11), OUT(100,12), F(11), G(10,11), BETA(11),
     * MARK(11), DXDT(11)
C OBTAINING OF DEPENDENCE OF SOLUTION X(ALFA) OF EQUATION
C F ( X , ALFA ) = 0
C ON PARAMETER ALFA BY MODIFIED METHOD OF DIFFERENTIATION
C WITH RESPECT TO PARAMETER
C N - NUMBER OF UNKNOWNS X(I)
C X(1),...,X(N) - INITIAL VALUES OF X(I), AFTER RETURN FINAL VALUES
C                 OF X(I)
C X(N+1) - INITIAL VALUE OF PARAMETER ALFA, AFTER RETURN FINAL VALUE
C          OF ALFA
C XLOW(1),...,XLOW(N) - LOWER BOUNDS FOR X(I)
C XUPP(1),...,XUPP(N) - UPPER BOUNDS FOR X(I)
C XLOW(N+1),XUPP(N+1) - LOWER AND UPPER BOUNDS FOR ALFA
C    IF XLOW OR XUPP IS EXCEEDED, THEN END OF DERPAR
C    AND MAXOUT=-2 AFTER RETURN
C EPS - ACCURACY DESIRED IN NEWTON ITERATION FOR
C    SUM OF (W(I)*ABS(XNEW(I)-XOLD(I))), I=1,...,N+1
C W(1),...,W(N+1) - WEIGHTS USED IN TERMINATION CRITERION OF NEWTON
C                   PROCESS
C INITAL - IF (INITAL.NE.0) THEN SEVERAL STEPS IN NEWTON ITERATION
C         ARE MADE BEFORE COMPUTATION IN ORDER TO INCREASE
C         ACCURACY OF INITIAL POINT-
C     IF (INITAL.EQ.1.AND.EPS-ACCURACY IS NOT FULFILLED IN ITIN
C     ITERATIONS) THEN RETURN.   IF (INITAL.EQ.2) THEN ALWAYS RETURN
C     AFTER INITIAL NEWTON ITERATION, RESULTS ARE IN X.
C     IF (INITAL.EQ.3) THEN CONTINUE IN DERPAR AFTER INITIAL NEWTON.
C     IF (INITAL.EQ.0) THEN NO INITIAL NEWTON ITERATION IS USED.
C ITIN - MAXIMAL NUMBER OF INITIAL NEWTON ITERATIONS. IF
C       EPS-ACCURACY IS NOT FULFILLED IN ITIN ITERATIONS THEN
C       ITIN=-1 AFTER RETURN.
C HH - INTEGRATION STEP ALONG ARC LENGTH OF SOLUTION LOCUS
C HMAX(1),...,HMAX(N+1) - UPPER BOUNDS FOR INCREMENTS OF X(I) IN
C                        ONE INTEGRATION STEP (APPROXIMATION ONLY)
C PREF(1),...,PREF(N+1) - PREFERENCE NUMBERS (EXPLANATION SEE IN
C                         SUBR.GAUSE)
C NDIR(1),...,NDIR(N+1) - INITIAL CHANGE OF X(I) IS POSITIVE ALONG
C     SOLUTION LOCUS (CURVE) IF NDIR(I)=1 AND NEGATIVE IF
C     NDIR(I)=-1.
C E - CRITERION FOR TEST ON CLOSED CURVE, IF
C    (SUM OF (W(I)*ABS(X(I)-XINITIAL(I))),I=1,...,N+1).LE.E)
C    THEN CLOSED CURVE MAY BE EXPECTED.
C MXADMS - MAXIMAL ORDER OF ADAMS-BASHFORTH FORMULA,
C          1.LE.MXADMS.LE.4.
C NCORR - MAXIMAL NUMBER OF NEWTON CORRECTIONS AFTER PREDICTION
C    BY ADAMS-BASHFORTH METHOD.
C NCRAD - IF (NCRAD.NE.0) THEN ADDITIONAL NEWTON CORRECTION
C    WITHOUT NEW COMPUTATION OF JACOBIAN MATRIX IS USED.
C NOUT - AFTER RETURN NUMBER OF CALCULATED POINTS ON THE CURVE
C        X(ALFA), NOUT.LE.MAXOUT.
C OUT(J,1),OUT(J,2),...,OUT(J,N+1) - J-TH POINT X(1),...,X(N),ALFA
C    ON CURVE X(ALFA)
C OUT(J,N+2) - VALUE OF SQRT(SUM OF SQUARES OF F).
C    IF (OUT(J,N+2).LT.0.0) THEN  ABS(OUT(J,N+2)) CORRESPONDS TO X
C    AND ALFA NOT EXACTLY, BECAUSE ADDITIONAL NEWTON CORRECTION WAS
C    USED (NCRAD.NE.0).
C    VALUES F(I) ARE NOT COMPUTED FOR X AND ALFA PRINTED AND/OR
C    STORED AND THEREFORE LAST TIME COMPUTED VALUE OF SQRT(SUM OF
C    SQUARES OF F) IS ON OUR DISPOSAL ONLY.
C MAXOUT - MAXIMAL NUMBER OF CALCULATED POINTS ON CURVE X(ALFA).
C    IF MAXOUT AFTER RETURN EQUALS TO-
C          -1 - THEN CLOSED CURVE X(ALFA) MAY BE EXPECTED
C          -2 - THEN BOUND XLOW OR XUPP WAS EXCEEDED
C          -3 - THEN SINGULAR JACOBIAN MATRIX OCCURRED, ITS RANK IS
C               .LT. N.
C NPRNT - IF (NPRNT.EQ.3) THEN RESULTING POINTS ON CURVE X(ALFA)
C    ARE IN ARRAY OUT( , ) AFTER RETURN. IF (NPRNT.EQ.1) THEN THESE
C    POINTS ARE PRINTED ONLY.  IF (NPRNT.EQ.2) THEN THESE POINTS ARE
C    BOTH PRINTED AND IN ARRAY OUT.
C SUBROUTINE FCTN MUST BE PROGRAMMED IN FOLLOWING FORM -
C   SUBROUTINE FCTN(N,X,F,G)
C   DIMENSION X(11),F(10),G(10,11)
C   F(I)= FI (X(1),X(2),...,X(N),ALFA) FOR I=1,...,N
C   G(I,J)= D FI (X(1),...,X(N),ALFA)/ D X(J) FOR  I,J=1,...,N
C   G(I,N+1)= D FI (X(1),...,X(N),ALFA)/ D ALFA  FOR I=1,...,N
C   RETURN
C   END
      DATA INDIC, INDSP /1H*,1H /
C LW IS PRINTER DEVICE NUMBER
      N1 = N + 1
      LW = 3
      IF (INITAL) 10, 60, 10
C INITIAL NEWTON ITERATIONS
   10 DO 40 L=1,ITIN
        CALL FCTN(N, X, F, G)
        SQUAR = 0.0
        DO 20 I=1,N
          SQUAR = SQUAR + F(I)**2
   20   CONTINUE
        LL = L - 1
        SQUAR = SQRT(SQUAR)
        IF (NPRNT.NE.3) WRITE (LW,99999) LL, (X(I),I=1,N1), SQUAR
        CALL GAUSE(N, G, F, M, 10, 11, PREF, BETA, K)
        IF (M.EQ.0) GO TO 310
        P = 0.0
        DO 30 J=1,N1
          X(J) = X(J) - F(J)
          P = P + ABS(F(J))*W(J)
   30   CONTINUE
        IF (P.LE.EPS) GO TO 50
   40 CONTINUE
      WRITE (LW,99998) ITIN
      ITIN = -1
      IF (INITAL.EQ.1) RETURN
   50 IF (NPRNT.NE.3) WRITE (LW,99997) (X(I),I=1,N1)
      IF (INITAL.EQ.2) RETURN
C AFTER INITIAL NEWTON ITERATIONS
   60 IF (NPRNT.NE.3) WRITE (LW,99996)
      KOUT = 0
      NOUT = 0
      MADMS = 0
      NC = 1
      K1 = 0
   70 CALL FCTN(N, X, F, G)
      SQUAR = 0.0
      DO 80 I=1,N
        SQUAR = SQUAR + F(I)**2
   80 CONTINUE
      CALL GAUSE(N, G, F, M, 10, 11, PREF, BETA, K)
      IF (M.EQ.0) GO TO 310
      IF (K1.EQ.K) GO TO 90
C CHANGE OF INDEPENDENT VARIABLE (ITS INDEX = K NOW)
      MADMS = 0
      K1 = K
   90 SQUAR = SQRT(SQUAR)
      IF (NCRAD.EQ.1) SQUAR = -SQUAR
      P = 0.0
      DO 100 I=1,N1
        P = P + W(I)*ABS(F(I))
  100 CONTINUE
      IF (P.LE.EPS) GO TO 130
      IF (NC.GE.NCORR) GO TO 120
      DO 110 I=1,N1
        X(I) = X(I) - F(I)
  110 CONTINUE
C ONE ITERATION IN NEWTON METHOD
      NC = NC + 1
      GO TO 70
  120 IF (NCORR.EQ.0) GO TO 130
      WRITE (LW,99995) NCORR, P
  130 NC = 1
      IF (NCRAD.EQ.0) GO TO 150
C ADDITIONAL NEWTON CORRECTION
      DO 140 I=1,N1
        X(I) = X(I) - F(I)
  140 CONTINUE
  150 NOUT = NOUT + 1
      DO 160 I=1,N1
        MARK(I) = INDSP
  160 CONTINUE
      MARK(K) = INDIC
      IF (NPRNT.EQ.3) GO TO 170
      WRITE (LW,99994)
      WRITE (LW,99993) (X(I),MARK(I),I=1,N1), SQUAR
  170 IF (NPRNT.EQ.1) GO TO 200
      IF (NOUT.LE.100) GO TO 180
      WRITE (LW,99992)
      RETURN
  180 DO 190 I=1,N1
        OUT(NOUT,I) = X(I)
  190 CONTINUE
      OUT(NOUT,N+2) = SQUAR
      GO TO 210
  200 IF (NOUT.EQ.1) GO TO 180
  210 IF (NOUT.GE.MAXOUT) RETURN
      DO 220 I=1,N1
        IF (X(I).LT.XLOW(I) .OR. X(I).GT.XUPP(I)) GO TO 300
  220 CONTINUE
      IF (NOUT.LE.3) GO TO 240
      P = 0.0
      DO 230 I=1,N1
        P = P + W(I)*ABS(X(I)-OUT(1,I))
  230 CONTINUE
      IF (P.LE.E) GO TO 290
C CLOSED CURVE MAY BE EXPECTED
  240 DXK2 = 1.0
      DO 250 I=1,N1
        DXK2 = DXK2 + BETA(I)**2
  250 CONTINUE
      DXDT(K) = 1.0/SQRT(DXK2)*FLOAT(NDIR(K))
C DERIVATIVE OF INDEPENDENT VARIABLE X(K) WITH RESPECT TO ARC LENGTH
C OF SOLUTION LOCUS IS COMPUTED HERE
      H = HH
      DO 270 I=1,N1
        NDIR(I) = 1
        IF (I.EQ.K) GO TO 260
        DXDT(I) = BETA(I)*DXDT(K)
  260   IF (DXDT(I).LT.0.0) NDIR(I) = -1
        IF (H*ABS(DXDT(I)).LE.HMAX(I)) GO TO 270
        MADMS = 0
        H = HMAX(I)/ABS(DXDT(I))
  270 CONTINUE
      IF (NOUT.LE.KOUT+3) GO TO 280
      IF (H*ABS(DXDT(K)).LE.0.8*ABS(X(K)-OUT(1,K))) GO TO 280
      IF ((OUT(1,K)-X(K))*FLOAT(NDIR(K)).LE.0.0) GO TO 280
      MADMS = 0
      IF (H*ABS(DXDT(K)).LE.ABS(X(K)-OUT(1,K))) GO TO 280
      H = ABS(X(K)-OUT(1,K))/ABS(DXDT(K))
      KOUT = NOUT
  280 CALL ADAMS(N, DXDT, MADMS, H, X, MXADMS)
      GO TO 70
  290 WRITE (LW,99991)
      MAXOUT = -1
      RETURN
  300 MAXOUT = -2
      RETURN
  310 WRITE (LW,99990) (X(I),I=1,N1)
      MAXOUT = -3
      RETURN
99999 FORMAT (3X, 7H DERPAR, I3, 25H.INITIAL NEWTON ITERATION/22X,
     * 11HX,ALFA,SQF=, 5F15.7/(33X, 5F15.7))
99998 FORMAT (/16H DERPAR OVERFLOW, I5, 2X, 18HINITIAL ITERATIONS/)
99997 FORMAT (3X, 39H DERPAR AFTER INITIAL NEWTON ITERATIONS/22X,
     * 7HX,ALFA=, 4X, 5F15.7/(33X, 5F15.7))
99996 FORMAT (/6X, 45H DERPAR RESULTS (VARIABLE CHOSEN AS INDEPENDE,
     * 18HNT IS MARKED BY *)/)
99995 FORMAT (/37H DERPAR NUMBER OF NEWTON CORRECTIONS=, I5,
     * 31H  IS NOT SUFFICIENT,ERROR OF X=, F15.7/)
99994 FORMAT (6X, 27H DERPAR RESULTS X,ALFA,SQF=)
99993 FORMAT (33X, 5(F15.7, A1))
99992 FORMAT (/28H DERPAR OUTPUT ARRAY IS FULL//)
99991 FORMAT (/36H DERPAR CLOSED CURVE MAY BE EXPECTED/)
99990 FORMAT (/48H DERPAR SINGULAR JACOBIAN MATRIX FOR X AND ALFA=,
     * 5F12.6/(48X, 5F12.6))
      END
      SUBROUTINE ADAMS(N, D, MADMS, H, X, MXADMS)                       ADA   10
      DIMENSION DER(4,11), X(11), D(11)
C ADAMS-BASHFORTH METHODS
      N1 = N + 1
      DO 20 I=1,3
        DO 10 J=1,N1
          DER(I+1,J) = DER(I,J)
   10   CONTINUE
   20 CONTINUE
      MADMS = MADMS + 1
      IF (MADMS.GT.MXADMS) MADMS = MXADMS
      IF (MADMS.GT.4) MADMS = 4
      DO 70 I=1,N1
        DER(1,I) = D(I)
        GO TO (30, 40, 50, 60), MADMS
   30   X(I) = X(I) + H*DER(1,I)
        GO TO 70
   40   X(I) = X(I) + 0.5*H*(3.0*DER(1,I)-DER(2,I))
        GO TO 70
   50   X(I) = X(I) + H*(23.0*DER(1,I)-16.0*DER(2,I)+5.0*DER(3,I))/
     *   12.0
        GO TO 70
   60   X(I) = X(I) + H*(55.0*DER(1,I)-59.0*DER(2,I)+37.0*DER(3,I)
     *   -9.0*DER(4,I))/24.0
   70 CONTINUE
      RETURN
      END
      SUBROUTINE GAUSE(N, A, B, M, NN, MM, PREF, BETA, K)               GAU   10
      DIMENSION A(NN,MM), B(MM), PREF(MM), BETA(MM), Y(11), X(11),
     * IRR(11), IRK(11)
C SOLUTION OF N LINEAR EQUATIONS FOR N+1 UNKNOWNS
C BASED ON GAUSSIAN ELIMINATION WITH PIVOTING.
C N - NUMBER OF EQUATIONS
C A - N X (N+1) MATRIX OF SYSTEM
C B - RIGHT-HAND SIDES
C M - IF (M.EQ.0) AFTER RETURN THEN RANK(A).LT.N
C PREF(I) - PREFERENCE NUMBER FOR X(I) TO BE INDEPENDENT VARIABLE,
C 0.0.LE.PREF(I).LE.1.0, THE LOWER IS PREF(I) THE HIGHER IS
C PREFERENCE OF X(I).
C BETA(I) - COEFFICIENTS IN EXPLICIT DEPENDENCES OBTAINED IN FORM -
C          X(I)=B(I)+BETA(I)*X(K), I.NE.K.
C K - RESULTING INDEX OF INDEPENDENT VARIABLE
      N1 = N + 1
      ID = 1
      M = 1
      DO 10 I=1,N1
        IRK(I) = 0
        IRR(I) = 0
   10 CONTINUE
   20 IR = 1
      IS = 1
      AMAX = 0.0
      DO 60 I=1,N
        IF (IRR(I)) 60, 30, 60
   30   DO 50 J=1,N1
          P = PREF(J)*ABS(A(I,J))
          IF (P-AMAX) 50, 50, 40
   40     IR = I
          IS = J
          AMAX = P
   50   CONTINUE
   60 CONTINUE
      IF (AMAX.NE.0.0) GO TO 70
C THIS CONDITION FOR SINGULARITY OF MATRIX MUST BE SPECIFIED
C MORE EXACTLY WITH RESPECT TO COMPUTER ACTUALLY USED
      M = 0
      GO TO 150
   70 IRR(IR) = IS
      DO 90 I=1,N
        IF (I.EQ.IR .OR. A(I,IS).EQ.0.0) GO TO 90
        P = A(I,IS)/A(IR,IS)
        DO 80 J=1,N1
          A(I,J) = A(I,J) - P*A(IR,J)
   80   CONTINUE
        A(I,IS) = 0.0
        B(I) = B(I) - P*B(IR)
   90 CONTINUE
      ID = ID + 1
      IF (ID.LE.N) GO TO 20
      DO 100 I=1,N
        IR = IRR(I)
        X(IR) = B(I)/A(I,IR)
        IRK(IR) = 1
  100 CONTINUE
      DO 110 K=1,N1
        IF (IRK(K).EQ.0) GO TO 120
  110 CONTINUE
  120 DO 130 I=1,N
        IR = IRR(I)
        Y(IR) = -A(I,K)/A(I,IR)
  130 CONTINUE
      DO 140 I=1,N1
        B(I) = X(I)
        BETA(I) = Y(I)
  140 CONTINUE
      B(K) = 0.0
      BETA(K) = 0.0
  150 RETURN
      END

.