p =.
\[Lambda] = 1;
k0 = (2 \[Pi])/\[Lambda];
U0 = HeavisidePi[x/p];
fU = InverseFourierTransform[U0, x, k];

p0 = Table[
   krange = 60 k0;
   \[Delta]k = 0.01 k0;
   zmax = 50;
   dU = Table[
     dfU = Table[fU E^(I Sqrt[k0^2 - k^2] z), {k, -krange, krange, \[Delta]k}];
     dim = Dimensions[dfU][[1]];
     RotateRight[Fourier[dfU], {Ceiling[dim/2]}]
     , {z, 0.001, zmax, \[Lambda]}];
   ListDensityPlot[Abs[dU[[All, 1 ;; -1 ;; 100]]], PlotRange -> All, ColorFunction -> "SunsetColors", DataRange -> {{-(k0/(2 \[Delta]k)), k0/(2 \[Delta]k)}, {0, zmax}}, LabelStyle -> {Black, Bold}, FrameLabel -> {"x/\[Lambda]", "z/\[Lambda]"}, PlotLabel -> "|E|"]
   , {p, 1, 50, 1}];
ListAnimate[p0]