program main
!=======================================================================
! Chapter 10, code 4: eigenvalues by Numerov shooting and bisection.
! Solve y'' = g(x,lambda)*y, with y(a) = y(b) = 0.
! Scan lambda for endpoint sign changes, then refine each bracket.
!
! Example: a uniform vibrating string, y'' = -lambda*y on [0,1].
! Exact eigenvalues: lambda_m = (m*pi/(b-a))**2, m = 1,2,... .
! Normalize each mode so that the integral of y**2 over [a,b] is 1.
!
! To solve another eigenvalue problem:
!   1. Edit the interval, grid, number of modes, and search controls.
!   2. Edit coefficient(x,lambda) below the driver.
!   3. Disable check_exact or replace the analytic comparisons below.
! This version assumes homogeneous fixed-value boundaries at BOTH ends
! and an equation without a first-derivative term.
!=======================================================================
implicit none
integer, parameter :: n = 200, modes = 3, data_unit = 20
integer i, m, iterations(modes)
double precision a, b, h, lambda_min, lambda_max, scan_step, tolerance
double precision x(0:n), values(modes), shapes(0:n,modes), residuals(modes)
double precision exact, exact_shape, errors(modes), pi, coefficient
logical check_exact
external coefficient

!-------------------- Problem and search: edit here ----------------------
a = 0.0d0
b = 1.0d0
lambda_min = 0.0d0
lambda_max = 100.0d0
scan_step = 1.0d0
tolerance = 1.0d-10             ! Maximum width of the final lambda bracket
check_exact = .true.
! Change n and modes in the declaration above.
! Increase lambda_max if the range contains too few modes.
! Make scan_step small enough to separate consecutive roots.
!------------------------------------------------------------------------
if (n < 2 .or. modes < 1 .or. b <= a) stop 'Use n >= 2, modes >= 1, and b > a.'
if (lambda_max <= lambda_min .or. scan_step <= 0.0d0 .or. tolerance <= 0.0d0) &
    stop 'Use an increasing search interval, positive scan step, and positive tolerance.'
h = (b-a)/dble(n)
do i = 0,n
    x(i) = a+dble(i)*h
end do
x(n) = b
pi = acos(-1.0d0)

write(*,'(a)') 'Eigenvalues: Numerov shooting and bisection'
write(*,'(a,i8)') 'Grid intervals:   ', n
write(*,'(a,f12.6)') 'Grid spacing:     ', h
write(*,'(a,es12.4)') 'Bracket tolerance:', tolerance
call find_modes(coefficient,n,x,modes,lambda_min,lambda_max,scan_step, &
                tolerance,values,shapes,residuals,iterations)

write(*,'(a)') ''
if (check_exact) then
    write(*,'(a6,2(1x,a18),1x,a12,1x,a10)') &
        'mode','numerical lambda','exact lambda','abs error','bisections'
else
    write(*,'(a6,1x,a18,1x,a10)') 'mode','numerical lambda','bisections'
end if
do m = 1,modes
    if (check_exact) then
        exact = (dble(m)*pi/(b-a))**2
        write(*,'(i6,2(1x,f18.10),1x,es12.4,1x,i10)') &
            m,values(m),exact,abs(values(m)-exact),iterations(m)
    else
        write(*,'(i6,1x,f18.10,1x,i10)') m,values(m),iterations(m)
    end if
end do

! One file for all eigenfunctions. Columns: x, mode 1, mode 2, ... .
! Comment lines record the eigenvalues and normalization convention.
open(unit=data_unit,file='ch10_numerov_eigenvalues.dat', &
     status='replace',action='write')
write(data_unit,'(a)') '# Numerov eigenfunctions: x y_1 y_2 ...'
write(data_unit,'(a)') '# Each mode has integral(y**2 dx) = 1 (trapezoidal rule).'
do m = 1,modes
    write(data_unit,'(a,i6,a,es24.16)') '# mode ',m,' lambda = ',values(m)
end do
write(*,'(a)') ''
write(*,'(a)') 'Eigenfunctions: columns are x, then the normalized modes.'
errors = 0.0d0
do i = 0,n
    write(*,'(f10.6)',advance='no') x(i)
    do m = 1,modes
        write(*,'(1x,es16.8)',advance='no') shapes(i,m)
        if (check_exact) then
            exact_shape = sqrt(2.0d0/(b-a))*sin(dble(m)*pi*(x(i)-a)/(b-a))
            errors(m) = max(errors(m),abs(shapes(i,m)-exact_shape))
        end if
    end do
    write(*,'(a)') ''
    write(data_unit,'(1000(1x,es24.16))') x(i),shapes(i,:)
end do
close(data_unit)
write(*,'(a)') ''
write(*,'(a6,1x,a18)') 'mode','endpoint residual'
do m = 1,modes
    ! These are endpoint values AFTER normalization. Do not force them to zero.
    write(*,'(i6,1x,es18.8)') m,residuals(m)
end do
if (check_exact) then
    write(*,'(a6,1x,a20)') 'mode','largest shape error'
    do m = 1,modes
        write(*,'(i6,1x,es18.8)') m,errors(m)
    end do
end if
write(*,'(a)') 'Data file: ch10_numerov_eigenvalues.dat'
end program main


double precision function coefficient(x,lambda)
!=======================================================================
! Supply g(x,lambda) in y'' = g(x,lambda)*y.
! Uniform string: g = -lambda.
! A dimensionless string with variable density rho(x) has g = -lambda*rho(x).
! Numerov requires this linear homogeneous form with no y' term.
!=======================================================================
implicit none
double precision x, lambda
coefficient = -lambda
end function coefficient


subroutine find_modes(fcn,n,x,modes,lambda_min,lambda_max,scan_step, &
                      tolerance,values,shapes,residuals,iterations)
!=======================================================================
! The shooting mismatch is y(b;lambda). Eigenvalues make it zero.
! Scan successive intervals for a sign change, then call bisection.
! A sign scan can miss two or more roots inside one interval; reduce
! scan_step and repeat when adapting this code to a different equation.
! Analytic mode numbers assume the scan finds the lowest modes in order.
!=======================================================================
implicit none
integer n, modes, found, count, iterations(modes)
double precision fcn, x(0:n), lambda_min, lambda_max, scan_step, tolerance
double precision values(modes), shapes(0:n,modes), residuals(modes)
double precision left, right, f_left, f_right, root, trial(0:n), mismatch
logical bracketed
external fcn
found = 0
left = lambda_min
call numerov(fcn,n,x,left,trial,f_left)
! Handle an eigenvalue exactly at the initial scan endpoint once.
if (f_left == 0.0d0) then
    found = 1
    values(found) = left
    call normalize_mode(n,x,trial)
    shapes(:,found) = trial
    residuals(found) = trial(n)
    iterations(found) = 0
end if

do while (found < modes .and. left < lambda_max)
    right = min(left+scan_step,lambda_max)
    if (right <= left) stop 'Scan step cannot advance lambda.'
    call numerov(fcn,n,x,right,trial,f_right)
    bracketed = (f_left < 0.0d0 .and. f_right > 0.0d0) .or. &
                (f_left > 0.0d0 .and. f_right < 0.0d0)
    if (f_right == 0.0d0 .or. bracketed) then
        if (f_right == 0.0d0) then
            root = right
            count = 0
        else
            call bisect(fcn,n,x,left,right,f_left,tolerance,root,count)
        end if
        found = found+1
        values(found) = root
        iterations(found) = count
        call numerov(fcn,n,x,root,trial,mismatch)
        call normalize_mode(n,x,trial)
        shapes(:,found) = trial
        residuals(found) = trial(n)
    end if
    left = right
    f_left = f_right
end do
if (found < modes) stop 'Too few modes found: enlarge the range or reduce scan_step.'
end subroutine find_modes


subroutine bisect(fcn,n,x,left,right,f_left,tolerance,root,iterations)
!=======================================================================
! Refine a sign-changing bracket. Stop on lambda-bracket width, not
! on mismatch magnitude, which depends on arbitrary shooting amplitude.
! Local copies preserve the scan endpoints in the calling subroutine.
!=======================================================================
implicit none
integer n, iterations
integer, parameter :: max_iterations = 100
double precision fcn, x(0:n), left, right, f_left, tolerance, root
double precision lo, hi, flo, mid, fmid, trial(0:n)
external fcn
lo = left
hi = right
flo = f_left
iterations = 0
do while (hi-lo > tolerance)
    if (iterations >= max_iterations) stop 'Bisection iteration limit reached.'
    mid = lo+(hi-lo)/2.0d0
    if (mid <= lo .or. mid >= hi) stop 'Lambda bracket cannot shrink; check tolerance.'
    call numerov(fcn,n,x,mid,trial,fmid)
    iterations = iterations+1
    if (fmid == 0.0d0) then
        root = mid
        return
    end if
    if ((flo < 0.0d0 .and. fmid > 0.0d0) .or. &
        (flo > 0.0d0 .and. fmid < 0.0d0)) then
        hi = mid
    else
        lo = mid
        flo = fmid
    end if
end do
root = lo+(hi-lo)/2.0d0
end subroutine bisect


subroutine numerov(fcn,n,x,lambda,y,mismatch)
!=======================================================================
! Numerov propagation for y'' = g(x,lambda)*y:
! (1-h**2*g(i+1)/12)*y(i+1) =
!   2*(1+5*h**2*g(i)/12)*y(i) - (1-h**2*g(i-1)/12)*y(i-1).
!
! y(0) = 0 and y(1) = h fix an arbitrary nonzero amplitude.
! No Taylor approximation of an initial derivative is needed here:
! the homogeneous linear recurrence makes the eigenvalue independent
! of that amplitude. Normalize the completed mode after root finding.
! This starting choice applies to the zero-value left boundary only.
!=======================================================================
implicit none
integer n, i
double precision fcn, x(0:n), lambda, y(0:n), mismatch
double precision g(0:n), h, c, denominator
external fcn
h = (x(n)-x(0))/dble(n)
c = h*h/12.0d0
do i = 0,n
    g(i) = fcn(x(i),lambda)
    if (.not. (abs(g(i)) <= huge(1.0d0))) stop 'Nonfinite Numerov coefficient.'
end do
y(0) = 0.0d0
y(1) = h
do i = 1,n-1
    denominator = 1.0d0-c*g(i+1)
    if (abs(denominator) <= epsilon(1.0d0)*max(1.0d0,abs(c*g(i+1)))) &
        stop 'Near-zero Numerov denominator: refine the grid.'
    y(i+1) = (2.0d0*(1.0d0+5.0d0*c*g(i))*y(i) &
               -(1.0d0-c*g(i-1))*y(i-1))/denominator
    if (.not. (abs(y(i+1)) <= huge(1.0d0))) &
        stop 'Nonfinite Numerov solution: check equation, interval, and grid.'
end do
mismatch = y(n)
end subroutine numerov


subroutine normalize_mode(n,x,y)
! Trapezoidal approximation to integral y**2 dx; keep the initial sign.
implicit none
integer n
double precision x(0:n), y(0:n), h, norm
h = (x(n)-x(0))/dble(n)
norm = sqrt(h*(0.5d0*y(0)**2+sum(y(1:n-1)**2)+0.5d0*y(n)**2))
if (norm <= 0.0d0 .or. .not. (norm <= huge(1.0d0))) &
    stop 'Cannot normalize the eigenfunction.'
y = y/norm
end subroutine normalize_mode