program numerov_bound
c	Dario Mitnik
c	17-May-2002
		
c....... Numerov method to find eigenvalues and eigenfunctions

c........ The example calculates the energies between an interval (Emin,Emax)
c	  and the wavefunctions for a unidimensional rectangular potential well, 
c	  having an Epot=V0. The center of the well is located at x=0
c	  and it has a half width of r0. 

c.........The Numerov algorithm is implemented forward from 
c	  -npts to imatch, and prints the wavefunction on numerleft.out
c	  The backward integration (npts to imatch) is printed on numeright.out	

c..........It is assumed that the distance (h) between points is 1, 
c	   for a physical problem it is necessary to scale h.


	implicit real*8(a-h,o-z)
	parameter (mxpts=1000)
	common /blockpot/Vpot(-mxpts:mxpts),V0,r0,npts,imatch
	
	eps = 1e-08		!! accuracy for eigenvalues

c....... opening output files
	open(unit=12,file='Vpot.dat',status='unknown')
	open(unit=13,file='numerleft.dat',status='unknown')
	open(unit=14,file='numeright.dat',status='unknown')
	
c....... fixed data, replace it with readings !!!
	mxstps = 100		!! maximum number of steps in iteration
	npts=1000
	V0=-0.001
	r0=500
	imatch = 500		!! point for matching the functions
	Emin=-0.001
	Emax=-0.00085
	print*,' depth of the potential well:  V=:',V0
	print*,' half width of the potential well:  r0=:',r0
	print*,' Energy Range:  Emin:',Emin,' to Emax=:',Emax	
	print*,' Point for matching the functions: imatch=:',imatch
		
c....... Calculate the Potential 
	do 20 i=-npts,npts
	    call Vpotcalc(i)
20	continue
	
c....... Bisection Algorithm for founding the energy
	do 50 i=1,mxstps
	  E = (Emin+Emax)/2.0			!! divide energy range
	  print*,' step :',i,'  Energy = :',E
	  rE = rmatchdiff(E)
	  rEmax = rmatchdiff(Emax)
	  if (rE*rEmax.gt.0) then
	      Emax = E
	  else
	      Emin = E
	  endif
	  if (abs(rE).lt.eps) go to 100
50	continue

	if (i.ge.mxstps) stop 'iteration fails '

100	print*,' Founded the eigenvalue E=:',E
	print*,' after ',i,' iterations'

	call printwave
		
	stop
	end
	
c
c **************************************************	  
c
	real*8 function rmatchdiff(E)
	implicit real*8(a-h,o-z)
	
	parameter (mxpts=1000)
	common /blockpot/Vpot(-mxpts:mxpts),V0,r0,npts,imatch
	common /blockwave/wave(-mxpts:mxpts)
	
c........ This function solves the differential equation
c         first from left to right (1 to 1500) --> r0left
c	  then from rigth to left (2000 to 1500) --> r0right
c	  and return the value rmatch = r0left - r0right

	
c....... calculate the wavefunction (forward) at point rmatch
	 call numerov(E,1)
	 r0right = wave(imatch)
c....... calculate the wavefunction (inwards) at point rmatch
	 call numerov(E,0)
	 r0left = wave(imatch)
	 
	rmatchdiff = r0left - r0right
	
	return
	end
		
c
c **************************************************	  
c
	subroutine numerov(E,iforward)
	implicit real*8(a-h,o-z)
	
c....... Numerov algorithm
c	If iforward=1: propagates from right to left
c	   iforward=0: inward propagation
	
	parameter (mxpts=1000)
	common /blockpot/Vpot(-mxpts:mxpts),V0,r0,npts,imatch
	common /blockwave/wave(-mxpts:mxpts)
	
c....... Initial Conditions:
	a0 = 0.0
	a1 = 1.e-5
c.......needs to be changed if the scale is changed
	h2 = 1.0			!! step interval  ( calculate it !!)
	
	if (iforward.eq.1) then
	  wave(-npts)=a0
	  wave(-npts+1)=a1
	  do 100 i=-npts+2,imatch
			f0 = E - Vpot(i)
			f1 = E - Vpot(i-1)
			f2 = E - Vpot(i-2)
      			aa= 2.0 - 10.0*h2*f1/12.0
      			bb= 1.0 + h2*f2/12.0
     			cc= 1.0 + h2*f0/12.0
     		        wave(i)=(aa*wave(i-1)-bb*wave(i-2))/cc
100	  continue
	else
	  wave(npts)=a0
	  wave(npts-1)=a1
	  do 200 i=npts-2,imatch,-1
			f0 = E - Vpot(i)
			f1 = E - Vpot(i+1)
			f2 = E - Vpot(i+2)
      			aa= 2.0 - 10.0*h2*f1/12.0
      			bb= 1.0 + h2*f2/12.0
     			cc= 1.0 + h2*f0/12.0
     		        wave(i)=(aa*wave(i+1)-bb*wave(i+2))/cc
200	  continue
	endif	
	
	return
	end
c
c **************************************************	  
c
	subroutine Vpotcalc(i)
	implicit real*8(a-h,o-z)
		
c....... Calculate the Potential
c	In this case is a Potential well

	parameter (mxpts=1000)
	common /blockpot/Vpot(-mxpts:mxpts),V0,r0,npts,imatch

	if (abs(i).gt.r0) then
c........       Outside the well
		Vpot(i) = 0
	else
c........       Inside the well
		Vpot(i) = V0
	endif
	
	write(12,*) i,' ',Vpot(i)
	
	return
	end

c
c **************************************************	  
c
	subroutine printwave
	implicit real*8(a-h,o-z)
	
c....... Print final wavefunction
	
	parameter (mxpts=1000)
	common /blockpot/Vpot(-mxpts:mxpts),V0,r0,npts,imatch
	common /blockwave/wave(-mxpts:mxpts)

	do 100 i=-npts,imatch
	  write(13,*) i,wave(i)
100	continue
	do 200 i=imatch+1,npts
	  write(14,*) i,wave(i)
200	continue

	return
	end