# - Solves the Schrodinger equation in one dimension, for a potential given      
#   by the function 'potential(x,xmax)', for x in the range (-xmax,xmax).        
# - The  integration of the wave functipn Psi(x) starts from boundary conditions 
#   Psi(-xmax)=0 and Psi(-xmax+h)=0.1, where h is the integration step                  
# - The boundary condition at x=xmax is taken to be Psi(xmax)=0.                 
# - Bisection is used to find the energy after course search for
#   two energies bracketing a solution).                                     

 function potential(v0::Float64,x::Float64,xmax::Float64)::Float64
    v::Float64=0
    if abs(x) < xmax/5
       v=v0
    end
    return v
 end 

 function normalize(n::Int64,h::Float64,psi)
    norm::Float64=psi[1]^2+psi[n]^2
    for i=2:n-3
      norm=norm+4*psi[i]^2+2*psi[i+1]^2
    end 
    norm=norm+4*psi[n-1]^2
    norm=1/(norm*h/3)^0.5
    psi.=psi.*norm
    return nothing
 end

 function numerov(v0::Float64,nx::Int64,xmax::Float64,ee::Float64,psi)

    h::Float64=xmax/nx
    h2::Float64=h^2
    h12::Float64=h2/12

    psi[1]=0
    psi[2]=0.1

    fn::Float64 = 2*(potential(v0,-xmax,xmax)-ee)        
    q0::Float64 = psi[1]*(1-h12*fn)
    fn = 2*(potential(v0,-xmax+h,xmax)-ee)        
    q1::Float64 = psi[2]*(1-h12*fn)

    for n=3:2*nx+1
       q2::Float64 = h2*fn*psi[n-1]+2*q1-q0
       fn = 2*(potential(v0,n*h-xmax,xmax)-ee)        
       psi[n] = q2/(1-h12*fn)
       q0=q1
       q1=q2
    end  
    normalize(2*nx+1,h,psi)

 end 

 function matchboundary(v0,nx,xmax,e1,de,eps,psi)

    numerov(v0,nx,xmax,e1,psi) 
    b1=psi[2*nx+1]  
    e2=e1
    for i=1:100
      e2=e2+de
      numerov(v0,nx,xmax,e2,psi)
      b2=psi[2*nx+1]
      if b1*b2 < 0 
         break
      end
      e1=e2
      b1=b2
    end
    println(e1,"  ",e2)

    while abs(e2-e1) > eps
       e3=(e1+e2)/2
       numerov(v0,nx,xmax,e3,psi)
       b3=psi[2*nx+1]
       if b3*b1 <= 0 
          e2=e3      
          b2=b3 
       else
          e1=e3
          b1=b3
       end
       println(e1,"  ",e2)
    end

    return (e1+e2)/2

 end

# Main program

 v0=5.0
 nx=200
 xmax=1.
 e0=0.
 de=0.1
 eps=0.0000001
 dx=xmax/nx

 psi = zeros(Float64,2*nx+1)

 e1=matchboundary(v0,nx,xmax,e0,de,eps,psi)
 xgrid=range(-xmax,xmax,length=2*nx+1)

 fi=open("psi.dat","w")
 for i=1:2*nx+1
    println(fi,(i-1)*dx-xmax,"  ",psi[i])
 end 

 using Plots
 display(plot(xgrid,psi))