SUBROUTINE FODE(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
     $                PAR,IPAR)
C
C SUBROUTINE  FODE  IS USED BY SUBROUTINE  STEPS  TO SPECIFY THE
C ORDINARY DIFFERENTIAL EQUATION  DY/DS = G(S,Y) , WHOSE SOLUTION
C IS THE ZERO CURVE OF THE HOMOTOPY MAP.  S = ARC LENGTH,
C YP = DY/DS, AND  Y(S) = (LAMBDA(S), X(S)) .
C
C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHOA, RHOJAC.
C
      DOUBLE PRECISION DDOT,DNRM2,S,SUM,YPNORM
      INTEGER I,IERR,IFLAG,IK,J,K,KP1,KPIV,LW,N,NFE,NP1
C
C *****  ARRAY DECLARATIONS.  *****
C
      DOUBLE PRECISION Y(N+1),YP(N+1),YPOLD(N+1),A(N),PAR(1)
      INTEGER IPAR(1)
C
C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
      DOUBLE PRECISION QR(N,N+1),ALPHA(N),TZ(N+1)
      INTEGER PIVOT(N+1)
C
C *****  END OF DIMENSIONAL INFORMATION.  *****
C
C
      NP1=N+1
      NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C
C COMPUTE THE JACOBIAN MATRIX, STORE IT IN QR.
C
      IF (IFLAG .EQ. -2) THEN
C
C  QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX )  .
C
        DO 30 K=1,NP1
          CALL RHOJAC(A,Y(1),Y(2),QR(1,K),K,PAR,IPAR)
30      CONTINUE
      ELSE
        CALL F(Y(2),TZ)
        IF (IFLAG .EQ. 0) THEN
C
C      QR = ( A - F(X), I - LAMBDA*DF(X) )  .
C
          DO 100 J=1,N
100       QR(J,1)=A(J)-TZ(J)
          DO 120 K=1,N
            CALL FJAC(Y(2),TZ,K)
            KP1=K+1
            DO 110 J=1,N
110         QR(J,KP1)=-Y(1)*TZ(J)
120       QR(K,KP1)=1.0+QR(K,KP1)
        ELSE
C
C   QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I ) .
C
140       DO 150 J=1,N
150       QR(J,1)=TZ(J)-Y(J+1)+A(J)
          DO 170 K=1,N
            CALL FJAC(Y(2),TZ,K)
            KP1=K+1
            DO 160 J=1,N
160         QR(J,KP1)=Y(1)*TZ(J)
170       QR(K,KP1)=1.0-Y(1)+QR(K,KP1)
        ENDIF
      ENDIF
C
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM.
210   CALL DCPOSE(N,N,QR,ALPHA,PIVOT,IERR,TZ,YP)
      IF (IERR .EQ. 0) GO TO 220
      IFLAG=4
      RETURN
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS.
220   TZ(NP1)=1.0
      DO 240 LW=1,N
        I=NP1-LW
        IK=I+1
        SUM=0.0
        DO 230 J=IK,NP1
230     SUM=SUM+QR(I,J)*TZ(J)
240   TZ(I)=-SUM/ALPHA(I)
      YPNORM=DNRM2(NP1,TZ,1)
      DO 260 K=1,NP1
        KPIV=PIVOT(K)
260   YP(KPIV)=TZ(K)/YPNORM
      IF (DDOT(NP1,YP,1,YPOLD,1) .GE. 0.0) GO TO 280
      DO 270 I=1,NP1
270   YP(I)=-YP(I)
C
C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN  YPOLD .
280   DO 290 I=1,NP1
290   YPOLD(I)=YP(I)
      RETURN
      END