shootf.f

      SUBROUTINE SHOOTF(NVAR,V1,V2,DELV1,DELV2,N1,N2,X1,X2,XF,EPS,H1,HMI
     *N,F,DV1,DV2)
      PARAMETER (NP=20)
      DIMENSION V1(N2),DELV1(N2),V2(N1),DELV2(N1),F(NVAR),DV1(N2),
     *    DV2(N1),Y(NP),F1(NP),F2(NP),DFDV(NP,NP),INDX(NP)
      EXTERNAL DERIVS,RKQC
      CALL LOAD1(X1,V1,Y)
      CALL ODEINT(Y,NVAR,X1,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
      CALL SCORE(XF,Y,F1)
      CALL LOAD2(X2,V2,Y)
      CALL ODEINT(Y,NVAR,X2,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
      CALL SCORE(XF,Y,F2)
      J=0
      DO 12 IV=1,N2
        J=J+1
        SAV=V1(IV)
        V1(IV)=V1(IV)+DELV1(IV)
        CALL LOAD1(X1,V1,Y)
        CALL ODEINT(Y,NVAR,X1,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
        CALL SCORE(XF,Y,F)
        DO 11 I=1,NVAR
          DFDV(I,J)=(F(I)-F1(I))/DELV1(IV)
11      CONTINUE
        V1(IV)=SAV
12    CONTINUE
      DO 14 IV=1,N1
        J=J+1
        SAV=V2(IV)
        V2(IV)=V2(IV)+DELV2(IV)
        CALL LOAD2(X2,V2,Y)
        CALL ODEINT(Y,NVAR,X2,XF,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
        CALL SCORE(XF,Y,F)
        DO 13 I=1,NVAR
          DFDV(I,J)=(F2(I)-F(I))/DELV2(IV)
13      CONTINUE
        V2(IV)=SAV
14    CONTINUE
      DO 15 I=1,NVAR
        F(I)=F1(I)-F2(I)
        F1(I)=-F(I)
15    CONTINUE
      CALL LUDCMP(DFDV,NVAR,NP,INDX,DET)
      CALL LUBKSB(DFDV,NVAR,NP,INDX,F1)
      J=0
      DO 16 IV=1,N2
        J=J+1
        V1(IV)=V1(IV)+F1(J)
        DV1(IV)=F1(J)
16    CONTINUE
      DO 17 IV=1,N1
        J=J+1
        V2(IV)=V2(IV)+F1(J)
        DV2(IV)=F1(J)
17    CONTINUE
      RETURN
      END