(*Root-finding code idea from https://mathematica.stackexchange.com/questions/16439/find-all-roots-of-an-interpolating-function-solution-to-a-differential-equation/16444#16444*)
Clear[findAllRoots]
SyntaxInformation[findAllRoots] = {"LocalVariables" -> {"Plot", {2, 2}}, "ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllRoots, HoldAll];
Options[findAllRoots] = Join[{"ShowPlot" -> False, PlotRange -> All}, FilterRules[Options[Plot], Except[PlotRange]]];
findAllRoots[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] := Module[{pl, p, x, localFunction, brackets}, localFunction = ReleaseHold[Hold[fn] /. HoldPattern[l] :> x];
  If[lmin != lmax, pl = Plot[localFunction, {x, lmin, lmax}, Evaluate@FilterRules[Join[{opts}, Options[findAllRoots]], Options[Plot]]];
   p = Cases[pl, Line[{x__}] :> x, Infinity];
   If[OptionValue["ShowPlot"], Print[Show[pl, PlotLabel -> "Finding roots for this function", ImageSize -> 200, BaseStyle -> {FontSize -> 8}]]], p = {}];
  brackets = Map[First, Select[(*This Split trick pretends that two points on the curve are "equal" if the function values have _opposite _ sign.Pairs of such sign-changes form the brackets for the subsequent FindRoot*) Split[p, Sign[Last[#2]] == -Sign[Last[#1]] &], Length[#1] == 2 &], {2}];
  x /. Apply[FindRoot[localFunction == 0, {x, ##1}] &, brackets, {1}] /. x -> {}]
(*8*)
\[Lambda] =.;
k0[\[Lambda]_] := (2 \[Pi])/\[Lambda];
n0 = 1; n1 = 1.1; \[Mu]0 = 1; c = 1; \[Omega][\[Lambda]_] := k0[\[Lambda]] c;
r0 = 2;
\[Beta][\[Lambda]_, kz_] := Sqrt[k0[\[Lambda]]^2 n1^2 - kz^2];
\[Sigma][\[Lambda]_, kz_] := Sqrt[kz^2 - k0[\[Lambda]]^2 n0^2];
dispersionTE[\[Lambda]_, kz_] := BesselJ[1, \[Beta][\[Lambda], kz] r0]/(\[Beta][\[Lambda], kz] BesselJ[0, \[Beta][\[Lambda], kz] r0]) + BesselK[1, \[Sigma][\[Lambda], kz] r0]/(\[Sigma][\[Lambda], kz] BesselK[0, \[Sigma][\[Lambda], kz] r0]);
Er[r_, root_] := 0; H\[Phi][r_, root_] := 0;
E\[Phi][r_, root_] := Piecewise[{{-I (\[Omega][\[Lambda]] \[Mu]0 )/\[Beta][\[Lambda], kz] BesselJ[1, \[Beta][\[Lambda], kz] r], r < r0}, {I (\[Omega][\[Lambda]] \[Mu]0)/\[Sigma][\[Lambda], kz] BesselJ[0, \[Beta][\[Lambda], kz] r0]/ BesselK[0, \[Sigma][\[Lambda], kz] r0] BesselK[1, \[Sigma][\[Lambda], kz] r], r >= r0}}] /. {kz -> root};
Hr[r_, root_] := Piecewise[{{I  BesselJ[1, \[Beta][\[Lambda], kz] r], r < r0}, {-I \[Beta][\[Lambda], kz]/\[Sigma][\[Lambda], kz] BesselJ[0, \[Beta][\[Lambda], kz] r0]/ BesselK[0, \[Sigma][\[Lambda], kz] r0] BesselK[1, \[Sigma][\[Lambda], kz] r], r >= r0}}] /. {kz -> root};
Hz[r_, root_] := Piecewise[{{BesselJ[0, \[Beta][\[Lambda], kz] r], r < r0}, {BesselJ[0, \[Beta][\[Lambda], kz] r0]/ BesselK[0, \[Sigma][\[Lambda], kz] r0] BesselK[0, \[Sigma][\[Lambda], kz] r], r >= r0}}] /. {kz -> root};

modes = Table[
   rootsTE = findAllRoots[dispersionTE[\[Lambda], kz], {kz, 0, 25}];
   rootsTE = Sort[rootsTE[[Flatten@Position[Evaluate[dispersionTE[\[Lambda], #] & /@ rootsTE], _?(Abs[#] < 1 &)]]] ];
   Column[{
     Style[StringForm["TE modes. \!\(\*FractionBox[\(\[Lambda]\), \(R\)]\)=``", NumberForm[N[\[Lambda]/r0], {3, 2}]], Black, Bold, FontSize -> 14],
     GraphicsRow[
      Table[DensityPlot[Norm[Hz[Sqrt[x^2 + y^2], rootsTE[[j]]]]^2 + Norm[Hr[Sqrt[x^2 + y^2], rootsTE[[j]]]]^2 + Norm[E\[Phi][Sqrt[x^2 + y^2], rootsTE[[j]]]]^2, {x, -1.5 r0, 1.5 r0}, {y, -1.5 r0, 1.5 r0}, PlotPoints -> 50, PlotRange -> All, ColorFunction -> "AvocadoColors", Frame -> False], {j, 1, Dimensions[rootsTE][[1]]}]  ]
     }, Alignment -> Center]
   , {\[Lambda], 0.3, 2.1, 0.03}];
ListAnimate[Reverse[modes], 2]