ClearAll["Global`*"]

[Theta] = 0 Degree;
[Omega]p = 5000;

[Xi]pzero[[Omega]_, k_] = [Omega]/(Sqrt[2] k Cos[[Theta]]);

[Lambda]p[k_] = k^2 Sin[[Theta]]^2;

[CapitalGamma]pzero[k_] =
  BesselI[0, [Lambda]p[k]] Exp[-[Lambda]p[k]];

Z[CapitalGamma]pzero[[Omega]_, k_] =
  I Sqrt[[Pi]]
    Exp[- ([Xi]pzero[[Omega], k])^2] (1 +
     Erf[I [Xi]pzero[[Omega], k]]);

pt[[Omega]_, k_] = ([Omega]p^2 [CapitalGamma]pzero[k])/
   k^2 (1 + ([Omega] Z[CapitalGamma]pzero[[Omega], k])/(
     Sqrt[2] k Cos[[Theta]]));

dispersion[[Omega]_, k_] := pt[[Omega], k]

dispersionrelation[
  k_] := (r = NSolve[ dispersion[[Omega], k] == 0, [Omega]];
  s = Evaluate[[Omega] /. r];
  Return[s])

Table[{[Theta], dispersionrelation[k]}, {k, 0.001, 1, 0.01}]