module global_vars
    integer :: N = 100
    real :: R = 15.0, V_0 = 4.0, a = 5.0
    save
end module

function V(x)
    use global_vars
    real :: V
    real, intent(in) :: x
    if (x < a) then
        V = -V_0
    else
        V = 0
    end if
    return
end function

function F(x, l)
    use global_vars
    real :: F
    real, intent(in) :: l, x
    F = V(x) + l * (l + 1) / (x * x)
    return
end function

program bcond
    use global_vars
    implicit none
    integer :: i, j, out
    real :: l, x, h, F, pi
    real, dimension(:), allocatable :: mBdiag, mBoff, mDdiag_inv, work, wr, wi
    real, dimension(:, :), allocatable :: mA, mV, mC, vl, vr, wavefun
    !memory allocation
    allocate(mBdiag(N))
    allocate(mBoff(N))
    allocate(mDdiag_inv(N))
    allocate(work(4 * N))
    allocate(wr(N))
    allocate(wi(N))
    allocate(mA(N, N))
    allocate(mV(N, N))
    allocate(mC(N, N))
    allocate(vl(N, N))
    allocate(vr(N, N))
    allocate(wavefun(N, N))
    !init variables
    out = 42
    l = 1.0
    h = R / N
    pi = 3.1415926535
    do i = 1, N
        x = h * i
        do j = 1, N
            mA(i, j) = 0.0
        end do
        mA(i, i) = 24.0 / h ** 2 + 10.0 * F(x, l)
        !mBdiag(i) = 10.0
        !mBoff(i) = 1.0
        if (i > 1) then
            mA(i, i - 1) = -12.0 / h ** 2 + F(x - h, l)
        end if
        if(i < N) then
            mA(i, i + 1) = -12.0 / h ** 2 + F(x + h, l)
        end if
    end do
    !print *, mA
    !calc V and D
    !D will overwrite mBdiag
    !call ssteqr('I', N, mBdiag, mBoff(1:N - 1), mV, N, work, i)
    do i = 1, N
        do j = 1, N
            mV(i, j) = sqrt(2.0 / (N + 1.0)) * sin(i * j * pi / (N + 1.0))
        end do
        mBdiag(i) = 10.0 + 2.0 * cos(i * pi / (N + 1.0))
    end do
    !calc D^(-1)
    do i = 1, N
        mDdiag_inv(i) = 1.0 / mBdiag(i)
    end do
    !calc matrix C = D^(-1) * V^(-1) * A * V
    do i = 1, N
        do j = 1, N
            mC(i, j) = mDdiag_inv(i) * mV(j, i)
        end do
    end do
    mC = matmul(mC, mA)
    mC = matmul(mC, mV)
    !calc eigenvalues and eigenvectors of C * y = lambda * y
    !wr - real parts of eigenvalues
    !vr - real part of eigenvectors
    call sgeev('N', 'V', N, mC, N, wr, wi, vl, N, vr, N, work, 4 * N, i)
    !f = V * y
    wavefun = matmul(mV, vr)
    !output
    print *, "E:"
    do i = 1, N
        if(wr(i) < 0.0 .and. wr(i) > -V_0) then
            print *, wr(i)
        end if
    end do
    open(unit = out, file = "output.txt")
    do i = 1, N
        if(wr(i) < 0.0 .and. wr(i) > -V_0) then
            do j = 1, N
                write(out, "(F15.5 F15.5)") h * j, wavefun(j, i)
            end do
            write(out, *) " "
        end if
    end do
    close(out)
    print *, "X"
    do i = 1, N
        if(wr(i) < 0.0 .and. wr(i) > -V_0) then
            print *, a * sqrt(V_0 - abs(wr(i)))
        end if
    end do
end program