program shoot_main
      implicit none
      integer n2,nvar,nvar1,kmax,kount,KMAXX,NMAX
      parameter (n2=1,nvar1=2,NMAX=50,KMAXX=200)
      real f(n2),V(n2),x1,x2,dxsav,xp(KMAXX),yp(NMAX,KMAXX)
      logical check
      common /caller/ x1,x2,nvar
      common /path/ kmax,kount,dxsav,xp,yp

      nvar = 2
      kmax = KMAXX
      dxsav = 0.01

      write(*,*) 'Please type in beginning and final x : x1 and x2'
      read (*,*) x1, x2
      write(*,*) 'What is the guess of the initial condition ?'
      read (*,*) V

      call newt(V,n2,check)

      write(*,*) 'n2 is ;', n2
      write(*,*) 'The final value of V is  :', V
      write(*,*) 'check of newt is', check
      write(*,*) 'koumt is', kount

      call pg_draw_path()
      end



      subroutine load(x1,V,y)
      real x1,V,y(2)
      y(1) = 0.0
      y(2) = V
      end

      subroutine score(x2,y,f)
      integer n2
      parameter (n2=1)
      real x2,y(2),f(n2)
      f(1) = y(1) - 0.0
      write(*,*) 'f value in subroutine score is :', f(1) 
      end


c     subroutine derivs(x,y,dydx)
c     real x,y(2),dydx(2)
c     dydx(1) = y(2)
c     dydx(2) = 2.0
c     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

c How to give BC for a cannon problem ?

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


      subroutine pg_draw_path()
      implicit none
      integer pgopen,i,j,nvar,nvar1,kmax,kount,KMAXX,NMAX
      parameter (nvar1=2,NMAX=50,KMAXX=200)
      real x1,x2,dxsav,xp(KMAXX),yp(NMAX,KMAXX)
c For pgplot
      real y_min(nvar1),y_max(nvar1),yy(KMAXX),y_extra
c 
      common /caller/ x1,x2,nvar
      common /path/ kmax,kount,dxsav,xp,yp

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