Perhaps you are looking for something like this?
using Roots
N = 10^3
G, d, ω = 9.8, 30.0, 2π .* range(0.05, 1, N) # don't collect the range!!
k0 = @. ω^2 / G / tanh(ω * sqrt(d / G))
dispersion(k, G, d, ω) = k * G * tanh(k * d) - ω^2 # no dots, just a scalar function
disproots(G, d, ω, k0) = find_zero(k->dispersion(k, G, d, ω), k0)
Benchmark:
1.7.2> @benchmark disproots.($G, $d, $ω, $k0)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 279.100 μs … 1.387 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 295.400 μs ┊ GC (median): 0.00%
Time (mean ± σ): 310.837 μs ± 51.768 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▆██▆▇▅▅▄▃▃▃▂▂▂▂▂▁▁▁ ▂
████████████████████████▇▇██▆▇▇▇█▆▇▇▇█▇▇▇▆▇▇▇▅▇▅▅▅▅▅▄▆▅▆▆▅▆▅ █
279 μs Histogram: log(frequency) by time 536 μs <
Memory estimate: 7.94 KiB, allocs estimate: 1.
Seems to be approx 50-100x faster than your code at N = 1000, and runtime scales linearly with N, while your code seems to scale like N^2.
It seems to solve a somewhat different problem, as @skleinbo says, you are solving a system of N equations, that just happen to be completely decoupled, so the results are equal up to ~10^-16.