You wrote your nroots! function in a classic Matlab-ish “vectorized” style where it assumes that z is an array, and allocates the result nr as an array too. If z is a scalar, then size(z) = () and nr is a zero-dimensional array (which is not the same thing as a scalar in Julia).
In this case, it doesn’t seem like there is any advantage in a vectorized style — even if z is an array, the different z’s don’t seem to share any computation. I would suggest just writing a function that works for scalar z and then broadcasting it as needed. If I understand your code correctly, this should be:
function nroots(p, z; eps=1e-05, maxiter=50)
for _ in 1:maxiter
z, zprev= z-p(z)/(dp(z) - prod(z)/(2*dp(z))), z
nr += abs(z-zprev)^2 > eps
(Note that neither this nor your original nroots! function should end with a !, by Julia convention, since neither of them mutate their arguments. Reassigning z is not mutation.)
Then you can broadcast this as nr = nroots.(p, z) over an array z and it will return an array of scalars.