program newton 

    implicit none
    real (kind=8), parameter :: tol = 1e-10_8
    real (kind=8) :: x0, x1, err
    integer :: iterations

    interface
        function f(x) result(y)
                real (kind=8) :: x, y
        end function f
        function df(x) result(y)
                real (kind=8) :: x, y
        end function df
    end interface

    x0 =   3.0 
    err = abs(x0)
    iterations = 0

    do while ( iterations < 100 .and. err > tol )
        print *, iterations, x0, f(x0), err
        iterations = iterations+1
        x1 = x0 - f(x0)/df(x0)
        err = abs(x1 - x0)
        x0 = x1
    enddo

    if ( iterations < 100 ) then
        print *, "Converged to ", x1, " in ", iterations, " iterations."
    else
        print *, "Did not converge in ",  iterations, " iterations."
    endif

end program newton

function f(x) result(y)
        real (kind=8) :: x, y
        y = (x**2 - 1.0_8)/3.0_8 
        return
end function f

function df(x) result(y)
        real (kind=8) :: x, y
        y = 2*x/3 
        return
end function df