        PROGRAM D9R13
C       Driver for routine MNEWT
        PARAMETER(NTRIAL=8,TOLX=1.0E-6,N=2,TOLF=1.0E-6,NP=15)
        DIMENSION X(NP),ALPHA(NP,NP),BETA(NP)
        L=1
20      CONTINUE
                   WRITE(*,'(/1X,A,I2)') 'Starting vector number',L
                   DO 11 I=1,N
                           WRITE(*,*)'X(',I,')=?'
                           READ(*,*)X(I)
11                 CONTINUE
                   DO 13 J=1,NTRIAL
                        CALL MNEWT(1,X,N,TOLX,TOLF)
                        CALL USRFUN(X,ALPHA,BETA)
                        WRITE(*,'(/1X,T5,A,T14,A,T29,A/)')'I','X(I)','F'
                        DO 12 I=1,N
                                WRITE(*,'(1X,I4,2E15.6)')I,X(I),-BETA(I)
12                      CONTINUE
                        WRITE(*,'(/1X,A)') 'press RETURN to continue...'
                        READ(*,*)
13                      CONTINUE
        L=L+1
        GOTO 20
        END
C
        SUBROUTINE USRFUN(X,ALPHA,BETA)
        PARAMETER(NP=15)
        DIMENSION ALPHA(NP,NP),BETA(NP),X(NP)
        ALPHA(1,1)=-2.0*X(1)
        ALPHA(1,2)=1.0
        ALPHA(2,1)=1./X(1)/X(1)
        ALPHA(2,2)=1.0
        BETA(1)=X(2)-X(1)*X(1)
        BETA(2)=X(2)-(X(1)+1.)/X(1)
        BETA(1)=-BETA(1)
        BETA(2)=-BETA(2)
        RETURN
        END
      SUBROUTINE MNEWT(NTRIAL,X,N,TOLX,TOLF)
      PARAMETER (NP=15)
      DIMENSION X(NP),ALPHA(NP,NP),BETA(NP),INDX(NP)
      DO 13  K=1,NTRIAL
        CALL USRFUN(X,ALPHA,BETA)
        ERRF=0.
        DO 11 I=1,N
          ERRF=ERRF+ABS(BETA(I))
11      CONTINUE
        IF(ERRF.LE.TOLF)RETURN
        CALL LUDCMP(ALPHA,N,NP,INDX,D)
        CALL LUBKSB(ALPHA,N,NP,INDX,BETA)
        ERRX=0.
        DO 12 I=1,N
          ERRX=ERRX+ABS(BETA(I))
          X(I)=X(I)+BETA(I)
12      CONTINUE
        IF(ERRX.LE.TOLX)RETURN
13    CONTINUE
      RETURN
      END
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END
