c-----------------------------------------------------------------------
c Solve the dynamics of predator-prey described by:
c   dxdt(1) =  alpha*x(1) - beta *x(1)*x(2)   ... prey
c   dxdt(2) = -gamma*x(2) + delta*x(1)*x(2)   ... predator
c Model constants:
c   alpha
c   beta
c   gamma
c   delta
c Variables:
c   t = time
c   n = number of equations = 2
c Program Note:
c   Link with RK4 of "Numerical Recipes"
c-----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------
      parameter (n=2)
      dimension x(n),dxdt(n)
      external DERIV
      common /c/alpha,beta,gamma,delta

c Model parameters -----------------------------------------------------
      write(*,*)'***Solve the predator-prey dynamic equations***'
      write(*,*)'  dxdt(1) =  alpha*x(1) - beta* x(1)*x(2) ... Prey'
      write(*,*)'  dxdt(2) = -gamma*x(2) + delta*x(1)*x(2) ... Predator'
      write(*,600)' alpha: '
        read(*,*)alpha
      write(*,600)' beta:  '
        read(*,*)beta
      write(*,600)' gamma: '
        read(*,*)gamma
      write(*,600)' delta: '
        read(*,*)delta

c Assign initial values for integration --------------------------------
      write(*,600)' x1(0): '
        read(*,*)x(1)
      write(*,600)' x2(0): '
        read(*,*)x(2)

c Integration limits ---------------------------------------------------
      step   =  0.1
      t      =  0.
      tmax   = 10.
      istart =  1
      iend   = nint((tmax-t)/step)

c Calculate the dynamic values by integration --------------------------
      write(*,*)'Results are written in PREDATOR.DAT'
      open(1,file='predator.dat',status='new')
      write(1,604)t,x(1),x(2)

      do 100 i=istart, iend
c       RK4 requires the evaluation of the initial dxdt ! -- Do not forget this step
          call deriv(t,x,dxdt)
c       Actual integration with routine from "Numerical Recipes"
          call RK4(x,dxdt,n,t,step,x,deriv)
c       End of integration step
          t=t+step
c       Save results to a file
          write(1,604)t,x(1),x(2)
100   continue

      close(1)

c Some formats ---------------------------------------------------------
600   format(a\)
604   format(1p3e13.5)
      end


c ----------------------------------------------------------------------
      subroutine deriv(t,x,dxdt)
c ----------------------------------------------------------------------
c Predator-prey dynamics
c Global variables ...
c   x(1) ... prey
c   x(2) ... predator
c   dxdt(n) ... time derivative of x variables
c   n ... number of state variables to integrate
c   t ... time
c ----------------------------------------------------------------------
      dimension x(1),dxdt(1)
      common /c/alpha,beta,gamma,delta

      dxdt(1) =  alpha*x(1) - beta *x(1)*x(2)
      dxdt(2) = -gamma*x(2) + delta*x(1)*x(2)

      return
      end


c ----------------------------------------------------------------------
c From "Numerical Recipes"  p533
c ----------------------------------------------------------------------
      SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
      PARAMETER (NMAX=10)
      DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
      HH=H*0.5
      H6=H/6.
      XH=X+HH
      DO 11 I=1,N
        YT(I)=Y(I)+HH*DYDX(I)
11    CONTINUE
      CALL DERIVS(XH,YT,DYT)
      DO 12 I=1,N
        YT(I)=Y(I)+HH*DYT(I)
12    CONTINUE
      CALL DERIVS(XH,YT,DYM)
      DO 13 I=1,N
        YT(I)=Y(I)+H*DYM(I)
        DYM(I)=DYT(I)+DYM(I)
13    CONTINUE
      CALL DERIVS(X+H,YT,DYT)
      DO 14 I=1,N
        YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14    CONTINUE
      RETURN
      END