program picard

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

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

    x0 = 0.5 
    err = x0
    iterations = 0

    do while ( iterations < 100 .and. err > tol )
        print *, iterations, x0, f(x0), err
        iterations = iterations+1
        x1 = f(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 picard

function f(x) result(y)
        real (kind=8) :: x, y
        y = sqrt(3.0_8*x + 1) - x 
        return
end function