C  PCPRB4.FOR  27 February 1991
C
      PROGRAM PCPRB4
C
C  Ermentrout-Aronson-Kopell problem.
C
C  The function F(X) is the right hand side of an ODE, and is of interest
C  in studying the equilibrium behavior of the system.
C
C  F(X) is of the form:
C
C  FX(1)=R*(1-GAMMA-R*R)+S*(-MU*SIN(PHI)+GAMMA*COS(PHI))
C  FX(2)=S*(1-GAMMA-S*S)+R*( MU*SIN(PHI)+GAMMA*COS(PHI))
C  FX(3)=Q*(R*R-S*S)-GAMMA*SIN(PHI)*(R*R+S*S)/(R*S)
C        +MU*COS(PHI)*(R*R-S*S)/(R*S)
C
C  GAMMA is the parameter of interest.
C  Q may be fixed at -2.0, MU at 1.0, for this problem.
C  R, S and PHI are the fundamental variables.
C
C  One family of solutions is:
C  (1, 1, 0, GAMMA) which becomes unstable as GAMMA increases from zero.
C  Two branches of asymmetric solutions diverge from this solution at one
C  point.
C
C  Another family is:
C  (1-2*GAMMA, 1-2*GAMMA, PI, GAMMA)
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=4)
      PARAMETER (LIW=NVAR+29)
      PARAMETER (LRW=29+(6+NVAR)*NVAR)
C
      EXTERNAL  DENSLV
      EXTERNAL  FXBARD
      EXTERNAL  FPBARD
      EXTERNAL  PITCON
C
      REAL      FPAR(2)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(1)
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   LOUNIT
      CHARACTER NAME*12
      REAL      PARMU
      REAL      PARQ
      REAL      RWORK(LRW)
      REAL      XR(NVAR)
C
      LOUNIT=6
C
C  Set parameters
C
      DO 4 I=1,LIW
        IWORK(I)=0
4       CONTINUE
      DO 5 I=1,LRW
        RWORK(I)=0.0
5       CONTINUE
      PARQ=-2.0
      PARMU=1.0
      FPAR(1)=PARQ
      FPAR(2)=PARMU
C
C  Set some entries of work arrays.
C
C  IWORK(1)=0 ; This is a startup
C  IWORK(2)=4 ; Use X(4) for initial parameter
C  IWORK(3)=0 ; Program may choose parameter index
C  IWORK(4)=0 ; Update jacobian every newton step
C  IWORK(5)=1 ; Seek target values for X(1)
C  IWORK(6)=0 ; Seek no limit points.
C  IWORK(7)=1 ; Control amount of output.
C  IWORK(8)=6 ; Output unit
C  IWORK(9)=0 ; Jacobian choice.
C
      IWORK(1)=0
      IWORK(2)=4
      IWORK(3)=0
      IWORK(4)=0
      IWORK(5)=1
      IWORK(6)=0
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=0
C
C  RWORK(1)=0.0001 ; Absolute error tolerance
C  RWORK(2)=0.0001 ; Relative error tolerance
C  RWORK(3)=0.01   ; Minimum stepsize
C  RWORK(4)=0.10   ; Maximum stepsize
C  RWORK(5)=0.02   ; Starting stepsize
C  RWORK(6)=1.0    ; Starting direction
C  RWORK(7)=1.0    ; Target value (Seek solution with X(3)=1)
C
      RWORK(1)=0.0001
      RWORK(2)=0.0001
      RWORK(3)=0.01
      RWORK(4)=0.10
      RWORK(5)=0.02
      RWORK(6)=1.0
      RWORK(7)=1.0
C
C  Set starting point:
C
      XR(1)=0.71055
      XR(2)=1.05053
      XR(3)=0.46070
      XR(4)=0.5
C
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PCPRB4'
      WRITE(LOUNIT,*)'PITCON sample program.'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Ermentrout-Aronson-Kopell function'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Step  Type of Point     '//
     *'X(1)         X(2)         X(3)          X(4)'
      WRITE(LOUNIT,*)' '
      I=0
      NAME='Start point  '
      WRITE(LOUNIT,'(1X,I3,2X,A12,2X,4G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
      DO 10 I=1,30
        CALL PITCON(FPBARD,FPAR,FXBARD,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C
        IF(IWORK(1).EQ.1)THEN
          NAME='Corrected    '
        ELSEIF(IWORK(1).EQ.2)THEN
          NAME='Continuation '
        ELSEIF(IWORK(1).EQ.3)THEN
          NAME='Target point '
        ELSEIF(IWORK(1).EQ.4)THEN
          NAME='Limit point  '
          ENDIF
C
        WRITE(LOUNIT,'(1X,I3,2X,A12,2X,4G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
        IF(IWORK(1).EQ.3)THEN
          WRITE(LOUNIT,*)'We have reached the desired point.'
          GO TO 20
          ENDIF
C
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR
          GO TO 20
          ENDIF
C
10      CONTINUE
C
20    CONTINUE
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Jacobians      ',IWORK(19)
      WRITE(LOUNIT,*)'Factorizations ',IWORK(20)
      WRITE(LOUNIT,*)'Solves         ',IWORK(21)
      WRITE(LOUNIT,*)'Functions      ',IWORK(22)
      STOP
      END
      SUBROUTINE FXBARD(NVAR,FPAR,IPAR,X,FX,IERROR)
C
C***********************************************************************
C
C  Evaluate the nonlinear function
C
      INTRINSIC COS
      INTRINSIC SIN
C
      INTEGER   NVAR
C
      REAL      COS
      REAL      FPAR(2)
      REAL      FX(NVAR)
      REAL      GAMMA
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      PARMU
      REAL      PARQ
      REAL      PHI
      REAL      R
      REAL      S
      REAL      SIN
      REAL      X(NVAR)
C
      PARQ=FPAR(1)
      PARMU=FPAR(2)
      R=X(1)
      S=X(2)
      PHI=X(3)
      GAMMA=X(4)
C
      FX(1)=R*(1.0-GAMMA-R*R) + S*(-PARMU*SIN(PHI)+GAMMA*COS(PHI))
      FX(2)=S*(1.0-GAMMA-S*S) + R*( PARMU*SIN(PHI)+GAMMA*COS(PHI))
      FX(3)=PARQ*(R*R-S*S)-GAMMA*SIN(PHI)*(R*R+S*S)/(R*S)
     *     +PARMU*COS(PHI)*(R*R-S*S)/(R*S)
      RETURN
      END
      SUBROUTINE FPBARD(NVAR,FPAR,IPAR,X,FPRIME,IERROR)
C
C***********************************************************************
C
C  Evaluate the jacobian of the nonlinear function.
C
      INTRINSIC COS
      INTRINSIC SIN
C
      INTEGER   NVAR
C
      REAL      COS
      REAL      FPAR(2)
      REAL      FPRIME(4,4)
      REAL      GAMMA
      INTEGER   IERROR
      INTEGER   IPAR(1)
      REAL      PARMU
      REAL      PARQ
      REAL      PHI
      REAL      R
      REAL      S
      REAL      SIN
      REAL      X(NVAR)
C
      PARQ=FPAR(1)
      PARMU=FPAR(2)
      R=X(1)
      S=X(2)
      PHI=X(3)
      GAMMA=X(4)
C
      FPRIME(1,1)=1.0-GAMMA-3.0*R*R
      FPRIME(1,2)=-PARMU*SIN(PHI)+GAMMA*COS(PHI)
      FPRIME(1,3)=S*(-PARMU*COS(PHI)-GAMMA*SIN(PHI))
      FPRIME(1,4)=-R+S*COS(PHI)
      FPRIME(2,1)=(PARMU*SIN(PHI)+GAMMA*COS(PHI))
      FPRIME(2,2)=1.0-GAMMA-3.0*S*S
      FPRIME(2,3)=R*( PARMU*COS(PHI)-GAMMA*SIN(PHI))
      FPRIME(2,4)=-S+R*COS(PHI)
      FPRIME(3,1)=2.0*PARQ*R - GAMMA*SIN(PHI)*(R*R-S*S)/(R*R*S)
     *     +PARMU*COS(PHI)*(R*R+S*S)/(R*R*S)
      FPRIME(3,2)=-2.0*PARQ*S - GAMMA*SIN(PHI)*(S*S-R*R)/(S*S*R)
     *     -PARMU*COS(PHI)*(R*R+S*S)/(S*S*R)
      FPRIME(3,3)=-GAMMA*COS(PHI)*(R*R+S*S)/(R*S)
     *     -PARMU*SIN(PHI)*(R*R-S*S)/(R*S)
      FPRIME(3,4)=-SIN(PHI)*(R*R+S*S)/(R*S)
      RETURN
      END