program Assignment_3

    implicit none
    real(4), allocatable, dimension(:) ::L_p, w, A_un, r_arr
    real(4), allocatable, dimension(:,:) :: A, z, L_b
    real(4) :: k, x, lambda, ierr, h, sum, step, nm
    integer :: N, i, j, o, matz, max, r, half
    matz =  1
    half = 15
    N = 10
    lambda = 1.
    o = N
    step = 1e-5
    allocate( L_p(0:N), A_un(0:2*half))
    max = 2 * half
    allocate(A(1:N,1:N), L_b(1:N,0:max), r_arr(0:max))
    k = 2.
    nm = N
    do j = 0,N
        do i = 0, N
            A(j,i) = 0.
        end do
    end do
    h = 1.

    do i = 0, max
        r_arr(i) = i * step
    end do

!

    do j = 1, N
        do r = 0, max
            if (j==1) then
                L_p(0) = 1
                L_b(j, r) = SQRT((lambda * GAMMA(j + 1. - 1.))/GAMMA(2. + j)) * r * lambda * L_p(j - 1) * EXP(-r * lambda / 2.)
            else if (j==2) then
                L_p(1) = 1. + k - lambda * r
                L_b(j, r) = SQRT((lambda * GAMMA(j + 1. - 1.))/GAMMA(2. + j)) * r * lambda * L_p(j - 1) * EXP(-r * lambda / 2.)
            else
                L_p(j-1) = ((2. * (j-1.) + 1. + k - lambda * r) * L_p(j-2) - (j + k - 1.) * L_p(j - 3)) / ( j +1. -1.)

                L_b(j, r) = SQRT((lambda * GAMMA(j + 1. - 1.))/GAMMA(2. + j)) * r * lambda * L_p(j - 1) * EXP(-r * lambda / 2.)
            end if
        end do
    end do

    do j = 1, N

        do i = 1, j
            if (i == 1 .OR. i == 2) then
                do r = 0, max
                    h = L_b(2,r) - L_b(1,r)
                    A_un(r) = L_b(j,r) * (-1./2. * (15./4. * L_b(i,r) - 77./6. * L_b(i + 1,r) + 107./6 * L_b(i + 2,r) &
                            - 13. * L_b(i + 3,r) + 61./12. * L_b(i + 4,r) - 5./6. * L_b(i + 5,r)) / h**2 - L_b(i,r)/r)

                end do
                call integrate(A_un, half, h, sum)
                A(j,i) = sum
            else if (i == N .OR. i == N-1) then
                do r = 0, max
                    h = L_b(2,r) - L_b(1,r)
                    A_un(r) = L_b(j,r) * (-1./2. * (15./4. * L_b(i,r) - 77./6. * L_b(i - 1,r) + 107./6 * L_b(i - 2,r) &
                            - 13. * L_b(i - 3,r) + 61./12. * L_b(i - 4,r) - 5./6. * L_b(i - 5,r)) / h**2 - L_b(i,r)/r)
                end do
                !write(*,*)A_un
                call integrate(A_un, half, h, sum)
                A(j,i) = sum
            else
                do r = 0, max
                    h = L_b(2,r) - L_b(1,r)
                    A_un= L_b(j,r) * (-1./2. * (-1./12. * L_b(i-2,r) + 4./3. * L_b(i-1,r) - 5./2. * L_b(i,r) &
                            + 4./3. * L_b(i + 1,r) - 1./12. * L_b(i + 2,r)) / h**2. - L_b(i,r)/r)
                end do
                call integrate(A_un, half, h, sum)
                A(j,i) = sum
                !write(*,*)"it got here"
            end if
        end do
    end do
    write(*,*)A(1,1:N)
    call rs(nm,o,a,w,matz,z,ierr)



end program Assignment_3

subroutine integrate(A_un, half, h, sum)
    implicit none
    real(4) :: A_un(0:2*half)
    real(4) :: sum, h
    integer :: half, i
    sum = 0
    do i = 0, half
        sum = sum + A_un(2*half-2) + 4. * A_un(2*half-1) + A_un(2*half)
    end do

    sum = h/3. * sum

end subroutine integrate