program odeint_loop

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

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,y_max,yy(KMAXX),x_temp(1001),y_temp(1001)
      real 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

      x_temp(1) = x1
      y_temp(1) = ystart(1)

      do 1000, i=1,1000
      x2tmp = x1 + i*(x2-x1)/1000
c     call odeint(ystart,nvar,x1,x2tmp,eps,h1,hmin,nok,nbad,derivs,rkqs)
      do j=1,nvar
        y_ini(j) = ystart(j)
      end do
      call odeint(y_ini,nvar,x1,x2tmp,eps,h1,hmin,nok,nbad,derivs,rkqs)
      write(*,*) 'Subroutine odeint completed.'
      write(*,*) 'nok is',nok,'        nbad is',nbad
      write(*,*) 'Final value of yi are:'
      write(*,*) (y_ini(j),j=1,nvar)
      y_temp(i+1)=y_ini(1)
      x_temp(i+1)=x2tmp
 1000 continue

c odeint does not pass out kount, we need to find out how many xp points.
c     do i=1,KMAXX
c       icount = i
c       if (xp(i).gt.x2) goto 101
c     end do
 101  write(*,*) 'Total number of points is',kount

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

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

c Start plotting
      if (pgopen('/xwin').le.0) stop
      call pgpap(6.0,0.75)
      y_extra = (y_max-y_min)/10.0
      call pgenv (x1,x2,y_min-y_extra,y_max+y_extra,0,0)
      call pgline(1001,x_temp,y_temp)
      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