      SUBROUTINE INITAL
C...
C...  THE DAVIDENKO DIFFERENTIAL EQUATION CAN BE WRITTEN IN MATRIX
C...  FORM AS
C...
C...     -  -       -  -       -
C...     J*DX/DT = -F, X(T0) = X0                                 (1)(2)
C...
C...  WHERE
C...
C...     -
C...     J              JACOBIAN MATRIX FOR THE NONLINEAR EQUATIONS
C...                    (NLE)
C...     -
C...     X              SOLUTION VECTOR FOR THE DAVIDENKO ODE
C...
C...     T              INDEPENDENT VARIABLE IN THE DAVIDENKO METHOD
C...
C...     T0             INITIAL VALUE OF T
C...
C...     -
C...     X0             INITIAL CONDITION VECTOR FOR EQUATION (1)
C...
C...                                    -                           -
C...  IN THE CODING IN SUBROUTINE DERV, J IS STORED IN ARRAY A, AND F IS
C...  STORED IN ARRAY B.  THE UNCOUPLING OF THE DERIVATIVES IN EQUATION
C...  (1) IS DONE IN SUBROUTINE DERV BY A CALL TO SUBROUTINES DECOMP AND
C...  SOLVE (LISTED AND DISCUSSED IN DETAIL IN FORSYTHE, G. E., M. A.
C...  MALCOLM AND C. B. MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPU-
C...  TATAIONS, CHAPTER 3, PRENTICE-HALL, ENGLEWOOD-CLIFFS, N.J., 1977).
C...  ANY OTHER ROUTINE(S) FOR LINEAR ALGEBRAIC EQUATIONS CAN BE USED IN
C...                                                    --   -
C...  PLACE OF DECOMP AND SOLVE (THE GENERAL PROBLEM IS AX = B).
C...
C...                                                      -   -
C...  THE PROGRAMMING OF EQUATION (1) IN TERMS OF VECTORS X, DX/DT AND
C...  -             -
C...  F, AND MATRIX J, IS INTENDED TO ILLUSTRATE THE GENERAL FEATURES OF
C...  DAVIDENKO*S METHOD, I.E., THESE VECTORS AND MATRIX CAN BE EXPANDED
C...  TO ACCOMMODATE A PROBLEM OF ANY SIZE.  HERE WE DIMENSION FOR
C...  TWO NLE PROVIDED BY L. T. BIEGLER:
C...
C...     X1**2 + X2**2 = 17
C...
C...     2*X1**(1/3) + X2**(1/2) = 4
C...
C...  DOUBLE PRECISION IS USED TO POSSIBLY ENHANCE THE PERFORMANCE
C...  OF DAVIDENKO*S METHOD WHEN THE JACOBIAN MATRIX IS POORLY
C...  CONDITIONED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(N=2)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      X(N)
     2      /F/   DXDT(N)
     3     /CV/    A(N,N),      B(N),   WORK(N),   IPVT(N),
     4               COND
C...
C...  INITIAL ESTIMATE OF THE SOLUTION
      IF(NORUN.EQ.1)THEN
         X(1)=1.0D0
         X(2)=1.0D0
      ELSE IF(NORUN.EQ.2)THEN
         X(1)=0.1D0
         X(2)=5.0D0
      END IF
C...
C...  INITIALIZE ALL OF THE MODEL CALCULATIONS
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(N=2)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      X(N)
     2      /F/   DXDT(N)
     3     /CV/    A(N,N),      B(N),   WORK(N),   IPVT(N),
     4               COND
C...
C...  JACOBIAN MATRIX
      A(1,1)=2.0D0*X(1)
      A(1,2)=2.0D0*X(2)
      A(2,1)=2.0D0*(1.0D0/3.0D0)*X(1)**(-2.0D0/3.0D0)
      A(2,2)=      (1.0D0/2.0D0)*X(2)**(-1.0D0/2.0D0)
C...
C...  COMPUTE THE RIGHT HAND SIDE VECTOR FOR THE DAVIDENKO DIFFERENTIAL
C...  EQUATION.  NOTE THAT THESE STATEMENTS CONSTITUTE A DEFINITION OF
C...  THE ELEMENTS OF B (OR VECTOR F IN EQUATION (1)).  THEY FOLLOW
C...  DIRECTLY FROM THE STATEMENT OF THE NONLINEAR EQUATIONS
      B(1)=-(      X(1)**2            +X(2)**2            -17.0D0)
      B(2)=-(2.0D0*X(1)**(1.0D0/3.0D0)+X(2)**(1.0D0/2.0D0)- 4.0D0)
C...
C...  SOLVE FOR THE VECTOR OF DERIVATIVES, DX/DT, IN EQUATION (1) (IN
C...  EFFECT, EQUATION (1) IS PREMULTIPLIED BY THE INVERSE OF JACOBIAN
C...  MATRIX J)
      CALL DECOMP(N,N,A,COND,IPVT,WORK)
C...
C...  CHECK THE CONDITION OF MATRIX A BEFORE COMPUTING THE DERIVATIVE
C...  VECTOR, DX/DT.  IF A IS NEAR SINGULAR, TERMINATE THE CURRENT RUN
      IF(COND.LT.1.0D+20)GO TO 3
      WRITE(NO,4)COND
4     FORMAT(1H ,//,
     1 30H THE CONDITION NUMBER OF A IS ,E10.3,33H SO THE CURRENT RUN IS
     1 TERMINATED)
      NSTOP=1
C...
C...  COMPUTE THE DERIVATIVE VECTOR DX/DT (INDIRECTLY VIA ARRAY B)
3     CALL SOLVE(N,N,A,B,IPVT)
C...
C...  TRANSFER ARRAY B TO ARRAY DXDT SO THE DERIVATIVE VECTOR, DX/DT,
C...  IS AVAILABLE IN COMMON/F/
      DO 2 I=1,N
      DXDT(I)=B(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(N=2)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      X(N)
     2      /F/   DXDT(N)
     3     /CV/    A(N,N),      B(N),   WORK(N),   IPVT(N),
     4               COND
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING THE EQUATION RESIDUALS AND
C...  CONDITION NUMBER
      DIMENSION F1R(11), F2R(11), CONP(11), TP(11)
C...
C...  MONITOR PROGRESS OF THE CALCULATION
      WRITE(*,*)' T = ',T
C...
C...  INITIALIZE OUTPUT
      IF(T.LT.0.001D0)THEN
C...
C...     PRINT HEADING FOR THE SOLUTION
         WRITE(NO,1)
1        FORMAT(9X,   'T',8X,'X1',    8X,'X2',/,
     +         10X,       4X,'DX1/DT',4X,'DX2/DT',/,
     +         10X,       6X,'RES1',  6X,'RES2',/,
     +         10X,       5X,'CONDITION NUMBER',/)
C...
C...     INITIALIZE COUNTER FOR PLOTTED SOLUTION
         IP=0
      END IF
C...
C...  PRINT THE NUMERICAL SOLUTION.  THE ELEMENTS OF B MUST FIRST BE
C...  RECOMPUTED SINCE B CONTAINS THE VECTOR DX/DT AFTER THE CALL TO
C...  SOLVE IN SUBROUTINE DERV
3     B(1)=-(      X(1)**2            +X(2)**2            -17.0D0)
      B(2)=-(2.0D0*X(1)**(1.0D0/3.0D0)+X(2)**(1.0D0/2.0D0)- 4.0D0)
C...
C...  PRINT THE SOLUTION, DERIVATIVES OF EQUATION (1) AND THE RESI-
C...  DUALS OF THE NONLINEAR EQUATIONS AT THE BEGINNING AND END OF
C...  EACH RUN
      IP=IP+1
      IF((IP.EQ.1).OR.(IP.EQ.11))THEN
         WRITE(NO,2)T,(X(I),I=1,N),(DXDT(I),I=1,N),(B(I),I=1,N),COND
2        FORMAT(F10.1,2F10.5,/,10X,2F10.5,/,10X,2F10.5,/,16X,E11.3,/)
      END IF
C...
C...  STORE THE RESIDUALS AND CONDITION NUMBER FOR PLOTTING
       F1R(IP)=DLOG10(DABS(B(1)))
       F2R(IP)=DLOG10(DABS(B(2)))
      CONP(IP)=DLOG10(DABS(COND))
      TP(IP)=T
C...
C...  PLOT THE LOG OF THE RESIDUALS OF THE NLE TO INVESTIGATE THE
C...  EXPONENTIAL CONVERGENCE OF DAVIDENKO*S ALGORITHM
      IF(IP.LT.11)RETURN
      CALL SPLOTS(1,IP,TP,F1R)
      WRITE(NO,11)
      CALL SPLOTS(1,IP,TP,F2R)
      WRITE(NO,12)
      CALL SPLOTS(1,IP,TP,CONP)
      WRITE(NO,13)
11    FORMAT(1H ,//,16X,'LOG10 OF RESIDUAL OF EQUATION 1 VS T')
12    FORMAT(1H ,//,16X,'LOG10 OF RESIDUAL OF EQUATION 2 VS T')
13    FORMAT(1H ,//,16X,'LOG10 OF CONDITION NUMBER VS T')
      RETURN
      END
