Function definitions and some tests

myRound[x_,n_]:=StringTake[ ToString[Round[x,0.1^n]]<>"000000000000000000",StringLength[ToString[Round[x,1]]]+1+n] (* no sprintf in Mathematica :( *)

norm[l_,m_] := 2^(l+6)Sqrt[(2 l+1)/4/Pi (l-m)!/(l+m)!] (* not used *)
psi[x_,y_,t_,a_,l_,m_] := Block[ {r,phi,CosTheta,tt,t2},
r=Sqrt[x^2+y^2]/a;
phi=ArcTan[x,y];
CosTheta=0;(* Equatoreal plane only, no Legendre polynomials used here*)
tt=t/a;
t2=tt^2;
r^l /(4 t2+(1+r^2-t2)^2)^((l+1)/2) *
Cos[m phi-(l+1)*ArcTan[1+r^2-t2,2*tt]]
]

L=10;
n1010=1/Abs[psi[Sqrt[L/(L+2)],0,0,1,L,L]]//N; (*maximum for l=m at Sqrt[L/(L+2)] *)
Plot[ n1010 psi[x,0,0,1,L,L],{x,0,5},PlotRange->Full]

tTest = 2.22;
DensityPlot[ n1010 psi[x,y,tTest,1,L,L],{x,-4,4},{y,-4,4},AspectRatio->Automatic,PlotRange->{-1,1},PlotPoints->200,ImageSize->640,Epilog->Text[Style["t="<>myRound[tTest,2],{White,20}],{3,3.5}],FrameStyle->20,
ColorFunction->(Blend[{{0,Cyan},{0.4,Blue},{0.5,Black},{0.6,Red},{1,Yellow}},#]&)]

(* So the DensityPlot is too rough, interpolation in ListDensityPlot is used instead. *)

(* frames = Table[
 DensityPlot[ n1010 psi[x,y,t,1,L,L],{x,-4,4},{y,-4,4},AspectRatio->Automatic,PlotRange->{-18,18},PlotPoints->50,ImageSize->480,Epilog->Text[Style["t="<>myRound[t,2],20],{3,3.5}],FrameStyle->20],
{t,-3.0,3.0,0.05}];
Export[ "psi1010.avi",frames ] 
*)

{L,n1010}
eps=1*^-10; (* to avoid sampling at the origin *)
data=Table[ n1010 psi[x+eps,y,2,1,L,L],{x,-4,4,0.1},{y,-4,4,0.1}];
({Max[#],Min[#]}&)[data//Flatten]

ListDensityPlot[Reverse[data],Mesh->None,InterpolationOrder->3,PlotRange->{-1,1},ColorFunction->(Blend[{{0,Cyan},{0.4,Blue},{0.5,Black},{0.6,Red},{1,Yellow}},#]&),ImageSize->480,FrameTicks->None,Frame->False,Epilog->Text[Style["t="<>myRound[2,2],{White,20}],{70,78}]]

ColorFunction -> Blend[ Black ]

(* Sample x-y plane at times from -3 .. 3 *)
dataT=Table[ n1010 psi[x+eps,y,t,1,L,L],{t,-3,3,0.05},{x,-4,4,0.1},{y,-4,4,0.1}];

timeT = Table[t,{t,-3,3,0.05}];

ProgressIndicator[Dynamic[k1],{0,Length[timeT]}] 
framesI = Table[ ListDensityPlot[Reverse[dataT[[k1]]] (*Reverse to correct y-axis orientation*),
Mesh->None,InterpolationOrder->3,PlotRange->{-1,1},
ColorFunction->(Blend[{{0,Cyan},{0.4,Blue},{0.5,Black},{0.6,Red},{1,Yellow}},#]&),
ImageSize->640,FrameTicks->None,Frame->False,
Epilog->Text[Style["t="<>myRound[timeT[[k1]],2],{White,20}],{70,78}]],{k1,1,Length[timeT]}
];

Export[ "psi1010bk3.avi",framesI]

Export[ "00a.jpeg",framesI[[33]]]

ColorFunction -> ThermometerColors

ProgressIndicator[Dynamic[k],{0,Length[timeT]}]  
framesJ=Table[ListDensityPlot[Reverse[dataT[[k]]] (*correct y-axis orientation*),
Mesh->None,InterpolationOrder->3,PlotRange->{-1,1},
ColorFunction->"ThermometerColors",
ImageSize->640,FrameTicks->None,Frame->False,
Epilog->Text[Style["t="<>myRound[timeT[[k]],2],{Black,20}],{70,78}]],{k,1,Length[timeT]}];

Export[ "L10M10b.avi",framesJ]
Export[ "L10M10b.jpeg",framesJ[[30]]]

framesJ[[30]]

Plot3D [ PlotTheme->{"Classic" }]

ProgressIndicator[Dynamic[t3],{-3,3}]  
frames3D = Table[Plot3D[ n1010 psi[x+eps,y,t3,1,L,L] , {x,-4,4},{y,-4,4},PlotTheme->{"Classic","ClassicLights"},PlotRange->{-1,1},MeshStyle->None,PlotPoints->150,ImageSize->640,TicksStyle->Directive[FontSize-> 20],Epilog->Text[Style["t="<>myRound[t3,2],{20}],{1,1}]],{t3,-3,3,0.05}];

Export[ "L10M10b-3D.avi",frames3D]
Export[ "L10M10-3D.jpeg",framesJ[[30]]]

frames3D[[30]]