ℰ[κ_,a_]=(1/(2 (a^2)))MathieuCharacteristicA[κ , ((a^2)/(π^2))10]

Plot[{ℰ[k,1],((k^2)/2)},{k,0,3.5}]

err=0.00000001;
Plot[{ℰ[0,a],ℰ[1-err,a],ℰ[1+err,a],ℰ[2-err,a],ℰ[2+err,a],ℰ[3-err,a]},{a,1,3},Axes-> False,Frame->True,FrameLabel->{"a","ℰ"}]

ℰ[1.5,2]

Plot[{MathieuCharacteristicA[κ ,0.1]},{κ,0,2.1},Axes->False,Frame->True]

ψ[a_,q_,x_]=MathieuC[a,q,x]

Plot[{((Abs[ψ[0.1,0.03,x]])^2)},{x,0,85},PlotRange->All]

Plot3D[ℰ[kx,1]+ℰ[ky,1],{kx,-3,3},{ky,-3,3},PlotPoints->100,Mesh->False]

Plot3D[((Abs[ψ[0,1,x]ψ[0,1,y]])^2),{x,0,50},{y,0,50},PlotPoints->50,Mesh->False,PlotRange->All,ViewPoint->{1.944, -2.504, 3.815}]

DensityPlot[((Abs[ψ[-0.2,1,x]ψ[-0.1,1,y]])^2),{x,0,10},{y,0,10},PlotPoints->100,Mesh->False,PlotRange->All]
DensityPlot[Cos[((2π x)/2)]+Cos[((2π y)/2)],{x,0,10},{y,0,10},PlotPoints->100,Mesh->False,PlotRange->All]