c
c This program is written by Prof. Ming-Hsien Lee for demonstrating the use
c of Numerical Recie subroutine ODEINT, RKQS, and RKCK to perform Adaptive
c Step-size Control Runge-Kutta algorithm for solving ordinary differential
c equations. Provided to year 2007 Numerical Methods class.
c
      program odeint_main

      implicit none
      integer nvar,i,nok,nbad
      parameter (nvar=2)
      real ystart(nvar),x1,x2,eps,h1,hmin

c For common block
      integer KMAXX,NMAX,kmax,kount
      parameter (KMAXX=200,NMAX=50)
      real xp(KMAXX),yp(NMAX,KMAXX),dxsav
      common /path/ kmax,kount,dxsav,xp,yp

c For pgplot
      integer pgopen,icount,j
      real y_min(nvar),y_max(nvar),yy(KMAXX),y_extra

      external rkqs,derivs

      kmax = KMAXX
      dxsav = 0.001

      write(*,*) 'What is the beginning and endding of x, x1 ; x2 ?'
      read (*,*) x1,x2
      write(*,*) 'Input inital vector of yi at x1:'
      read (*,*) (ystart(i),i=1,nvar)
      write(*,*) 'What is the accuracy eps and first setp size h1 ?'
      read (*,*) eps, h1 
      write(*,*) 'Minimum allowed step size (can be zero) hmin is:'
      read (*,*) hmin

      call odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs)
      write(*,*) 'Subroutine odeint completed.'
      write(*,*) 'Final value of y_i are:'
      write(*,*) (ystart(i),i=1,nvar)
      write(*,*) 'nok is',nok,'        nbad is',nbad
      write(*,*) 'The total number of points for plotting is :',kount

c Find Max and min values of yi(x), but here we will use only y1(x).
      do i=1,nvar
        y_max(i)= yp(i,1)
        y_min(i)= yp(i,1)
      end do
      do i=1,nvar
        do j=1,kount
          if (yp(i,j).gt.y_max(i)) y_max(i)=yp(i,j)
          if (yp(i,j).lt.y_min(i)) y_min(i)=yp(i,j)
        end do
      end do

c Copy y(1,1:nstep) to yy(1:nstep) for PGLINE plotting
      do i=1,kount
        yy(i) = yp(1,i)
      end do
c Intermidiate reult check
c     do i=1,kount
c       write(*,200) i,xp(i),i, yy(i)
c200    format(1x,'x(',i3,') is ',f5.2,' y(',i3,') is ',f5.2)
c     end do

c Start plotting
      y_extra = 0.1*(y_max(1)-y_min(1))
      if (pgopen('/xwin').le.0) stop
      call pgpap(6.0,0.75)
      call pgenv (x1,x2,y_min(1)-y_extra,y_max(1)+y_extra,0,0)
      call pgline(kount,xp,yy)
      call pgsci(2)
      call pgpt(kount,xp,yy,9)
      call pgclos

      end


      subroutine derivs(x,y,dydx)
      real x,y(2),dydx(2)
      dydx(1) = y(2)
      dydx(2) = -y(1)
      end

c     subroutine derivs(x,y,dydx)
c     real x,y,dydx
c     dydx = 2*x
c     end