solve_polynomial Function

pure function solve_polynomial(B, x0, error) result(x)

    real(wp),dimension(0:6),intent(in) :: B !! Polynomial `B = B(6) x^6 + B(5) x^5 + ... + B(0)`
    real(wp),intent(in) :: x0 !! Initial point
    real(wp),intent(in) :: error !! Maximum error
    real(wp) :: x !! root found after applying Newton-Raphson method to `B`
                  !! The function returns the value when the correction
                  !! is smaller than error.

    real(wp) :: f,fp,corr
    integer :: i, j !! counter

    integer,parameter :: maxiter = 100 !! maximum number of iterations

    x = x0
    do i = 1, maxiter
        f  = B(6)
        do j = 5, 0, -1
            if (j==5) then
                fp = f
            else
                fp = x*fp + f
            end if
            f  = x*f + B(j)
        end do
        if (fp==zero) exit ! singular point
        corr = f/fp
        x = x - corr
        if (abs(corr)<=error) exit
    end do

end function solve_polynomial