Program numerov
!-----------------------------------------------------------------
! April 2009
! 1D Schrodinger ODE y"(x)+2[E-V(x)]y(x)=0
! Numerical method: the Numerov's method for 2-nr oder ODE
! Root finding (e): the bisectional method
!-----------------------------------------------------------------
implicit none
integer, parameter :: n=101
double precision, dimension(-n:n) :: x, wl, wr
double precision e, a, b, emin, emax, energy, eps, fe, dx
double precision xmin, xmax
integer iflag, i
common/wf/ x, wl, wr
common/eqn/e, a, b, dx
external fe

open (unit=7,file="wfunc.dat")

eps = 1.0e-6
! Potential: well with a - width, b - depth
a =    1.0
b =   20.0
! domain for the wavefunction
xmax = 4.0
xmin = (-1.0)*xmax
x(n)  =  xmax
x(-n) =  xmin

! assumptions for two intial eigenvalues
emin = 2.0
emax = 4.0

! generate x(i) from x(-n), x(n) and n
dx = (x(n)-x(-n))/(2*n)
do i=-n,n,1
  x(i) = x(-n) + (i+n)*dx
!  write(*,100) i, x(i)
end do
100 format(i5,3f12.5)

!do while (emin <= 20.0)
! the bisectional search for a root for energy
  call Root_bs(fe,emin,emax,eps,energy,iflag)
  if(iflag.eq.1) then
    write(*,101) energy
  end if
!  emin = emin+2.0
!  emax = emax+2.0
!end do

101   Format(' Root is =',f12.5)
102   Format(' no root for Bisectional method')

write(7,200)
do i=-n,n,1
  write (7,201) x(i), wl(i), wr(i)
end do
200 format(5x,'x',11x,'y',11x,'z')
201	format(3(1pe12.4))
!102 format(' eigenvalue = ',1pe12.4)
stop
end

  function fe(energy)
!=================================================
integer, parameter :: n=101
double precision, dimension(-n:n) :: x, wl, wr
double precision fe
double precision e, energy, c1, c2, c3
double precision a, b, h
common/wf/ x, wl, wr
common/eqn/e, a, b, h
integer i, j, match

e=energy

! find a matching point (V(x)=energy)
do i=-n,n
  if (i >0 .and. g(x(i)) > 0.0) then
    match=i
	exit
  end if
end do
! for a step potential this point is x(match) = a
!do i=-n,n
!  if(x(i) > a) then
!    match = i
!	exit
!  end if
!end do
match = match
!write(*,200) match, x(match)
200 format (' match point = ',i5,2f12.4)

!define asympototic conditions (boudary values)
wl(-n)  = exp(-1.0*sqrt(2*abs(e))*abs(x(-n)))
wl(-n+1)= exp(-1.0*sqrt(2*abs(e))*abs(x(-n+1)))
wr(n)   = exp(-1.0*sqrt(2*abs(e))*abs(x(n)))
wr(n-1) = exp(-1.0*sqrt(2*abs(e))*abs(x(n-1)))

! integrate from the left (to the match point)
do i=-n+1,n-1,1
  c3 = 1.0+g(x(i+1))*h*h/12.0
  c2 = 2.0-(5.0/6.0)*h*h*g(x(i))
  c1 = 1.0+g(x(i-1))*h*h/12.0
  wl(i+1) = (c2*wl(i) - c1*wl(i-1))/c3
!  write(*,300) i+1, x(i+1), wl(i+1)
end do

! integrate from the right to the point im
do i=n-1,1-n,-1
  c3 = 1.0+g(x(i+1))*h*h/12.0
  c2 = 2.0-(5.0/6.0)*h*h*g(x(i))
  c1 = 1.0+g(x(i-1))*h*h/12.0
  wr(i-1) = (c2*wr(i) - c3*wr(i+1))/c1
!  write(*,300) i-1, x(i-1), wr(i-1)
end do

! rescale the right solution
do i=n-1,-n+1,-1
  wr(i) = wr(i)*(wl(match)/wr(match))
end do

fe = (wl(match+1)-wl(match-1))/(2.0*h*wl(match))
fe = fe - (wr(match+1)-wr(match-1))/(2.0*h*wr(match))
!write(*,103) energy, fe

100 format(i5,3f12.5)
300 format(i5,3e12.3)
103 format(f14.4,e14.4)

return
end

  Function g(x)
!--------------------------------------
! g(x) = 2*(E-V(x))
!--------------------------------------
  implicit none
  double precision g, x, a, b, e
  common/eqn/e, a, b
!  if(abs(x).le.a) then
!    g = 2.0*(e-b)
!	else
!	g = 2.0*e
!  end if
! potential from TP
  a = 1.0
  b = 4.0
  g = 0.5*a*a*b*(b-1.0)*(0.5-1.0/(cosh(a*x))**2)
  g = 2.0*(e-g)
  return
  end


        Subroutine root_bs(f,a,b,eps,Root,iflag)
!====================================================================
! Program to find a single root of an equation f(x)=0
! using the Bisectional method
! written by: Alex Godunov (October 2006)
!--------------------------------------------------------------------
! input ...
! f	- function which evaluates f(x) for any x in the interval a,b
! a	- left endpoint of initial interval
! b	- right endpoint of initial interval
! eps	- desired uncertainity of the root
!
! output ...
! Root  - root of equation f(x)=0
! iflag - indicator of success
!         1 - a single root found
!         0 - no solutions for the Bisectional method
!
! Limitations: a function f(x) must change sign between a and b
!
! Max number of allowed iterations is 100 (accuracy (b-a)/2**100)
!====================================================================
	Real*8 f,a,b,eps,Root
	Real*8 xl,x0,xr
	Integer*4 i, iter, iflag
	Parameter (iter=100)
!* check the bisection condition
	if(f(a)*f(b).gt.0.0) then
           iflag = 0
           return
	end if
	i=0
	xl = a
	xr = b
	do while (abs(xr-xl).gt.eps)
           i=i+1
           x0=(xr+xl)/2
           if(f(xl)*f(x0).le.0.0) then
              xr = x0
              else
              xl=x0
           end if
           if(i.ge.iter) exit
	end do
	Root=x0
	iflag=1
	return
    end