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)
dp= derivative(p)
ddp=derivative(dp)
prod=ddp*p
nr=0
for _ in 1:maxiter
z, zprev= z-p(z)/(dp(z) - prod(z)/(2*dp(z))), z
nr += abs(z-zprev)^2 > eps
end
return nr
end

(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.