      PROGRAM MAINPH
C
C     ------------------------------------------------------------
C
C     Testexample for PERHOM - Lorenz-Rayleigh-Benard-Problem
C
C*  Purpose           Testexample for code PERHOM
C*  Version           0.5 (Test Version)
C*  Revision          August 1984
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C
C     ------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION Z(3),Y(3,11),X(11),RW(6000),IW(100),ICALL(11),USCAL(3)
     1          ,IPAR(10),RH(2000),IH(100),FM(3,2)
      EXTERNAL PEDFXP,FCN
      COMMON /DATA/SIGMA,B
C
C  NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS
      N=3
C  NUMBER OF SHOOTING NODES
      M=11
C
      M1=M-1
C
C  AMOUNT OF REAL- AND INTEGER-WORKSPACE
      NRW=6000
      NIW=100
C
C  SHOOTING NODES
      X(1)=0.D0
      RED=1.D0/DBLE(M1)
      DO 1 I=2,M
1     X(I)=X(I-1)+RED
C
      SIGMA=1.6D1
      B=4.D0
C
C  INITIAL VALUES OF TRAJECTORY
      Y(1,1)=-0.06755D0
      Y(2,1)=-56.4816D0
      Y(3,1)=228.439D0
C
C  ESTIMATE OF PERIOD
      P=0.739342D0
C
C  VALUES OF HOMOTOPY PARAMETER AT START (PARS) AND AT END (PARE)
      PARS=242.000D0
      PARE=240.D0
C  INITIAL HOMOTOPY STEPSIZE
      SIG=0.005D0
C
C CALLING PEDFXP TO GENERATE INITIAL TRAJECTORY
      NRH=2000
      NIH=100
      IPAR(1)=0
      IPAR(2)=0
      IPAR(3)=0
      IPAR(4)=N
      IPAR(5)=0
      IPAR(6)=0
      T=0.D0
      H=1.D-3
      Z(1)=Y(1,1)
      Z(2)=Y(2,1)
      Z(3)=Y(3,1)
      TOL=1.D-3
      DO 22 I=2,M
      TEND=X(I)*P
      HMAX=TEND-T
      CALL PEDFXP (N,FCN,T,Z,TEND,PARS,TOL,HMAX,H,HS,USCAL,NRH,RH,
     1            NIH,IH,IPAR)
      Y(1,I)=Z(1)
      Y(2,I)=Z(2)
      Y(3,I)=Z(3)
22    CONTINUE
C
C  DESIRED RELATIVE ACCURACY FOR SOLUTION
      EPS=1.D-5
C  CLASSIFICATION OF RIGHT-HAND SIDE (AUTONOMOUS SYSTEM)
      ICALL(1)=1
C  PROBLEM IS NONLINEAR
      ICALL(2)=1
C  MAXIMUM PERMITTED GAUSS-NEWTON ITERATIONS PER HOMOTOPY STEP
      ICALL(3)=10
C  SOLVING VARIATIONAL EQUATION
      ICALL(4)=2
C  RANK-1 UPDATES ALLOWED
      ICALL(5)=1
C  ITERATIVE REFINEMENT IS ACTIVATED
      ICALL(6)=1
C  THE SYSTEM IS NONSTIFF
      ICALL(7)=0
C  PRINT PARAMETER
      ICALL(8)=2
C  FIRST TANGENT DIRECTION (TOWARDS PARE)
      ICALL(9)=-1
C  MAXIMUM PERMITTED HOMOTOPY STEPS
      ICALL(10)=20
C  SOLUTION DATA AT EACH STEP WILL NOT BE STORED
      ICALL(11)=0
C
C
      CALL PERHOM (PEDFXP,N,M,X,Y,P,PARS,PARE,SIG,EPS,FM,
     1                              NRW,RW,NIW,IW,ICALL)
C
C
      STOP
      END
C
      SUBROUTINE FCN(T,Z,PAR,DZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(3),DZ(3)
      COMMON /DATA/SIGMA,B
      DZ(1)=SIGMA*(Z(2)-Z(1))
      DZ(2)=Z(1)*(PAR-Z(3))-Z(2)
      DZ(3)=Z(1)*Z(2)-B*Z(3)
      RETURN
      END
C
      SUBROUTINE DFDY(T,Z,PAR,DF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(3),DF(3,3)
      COMMON /DATA/SIGMA,B
      DATA ZERO/0.D0/,ONE/1.D0/
      DF(1,1)=-SIGMA
      DF(1,2)=SIGMA
      DF(1,3)=ZERO
      DF(2,1)=PAR-Z(3)
      DF(2,2)=-ONE
      DF(2,3)=-Z(1)
      DF(3,1)=Z(2)
      DF(3,2)=Z(1)
      DF(3,3)=-B
      RETURN
      END
C
      SUBROUTINE DFDP(T,Y,PAR,DY)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION Y(3),DY(3)
      COMMON /DATA/SIGMA,B
      DY(1)=0.D0
      DY(2)=Y(1)
      DY(3)=0.D0
      RETURN
      END
