\[Lambda]0 = 1.; k0 = N[(2 \[Pi])/\[Lambda]0]; (*The wavelength in vacuum is set to 1, so all lengths are now in units of wavelengths*)
\[Delta] = \[Lambda]0/20; \[CapitalDelta] = 40*\[Lambda]0; (*Parameters for the grid*)

ReMapC[x_] := RGBColor[(2 x - 1) UnitStep[x - 0.5], 0, (1 - 2 x) UnitStep[0.5 - x]];
R = \[CapitalDelta]/3;
ren = Table[
   If[x^2 + y^2 <= R^2, 2/(1 + ((x^2 + y^2)/R^2)), 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
d = \[Lambda]0/2; (*typical scale of the absorbing layer*)
imn = Table[
   Chop[5 (E^-((x + \[CapitalDelta]/2)/d) + E^((x - \[CapitalDelta]/2)/d) + E^-((y + \[CapitalDelta]/2)/d) + E^((y - \[CapitalDelta]/2)/d))], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*Imaginary part of the refractive index (used to emulate absorbing boundaries)*)
dim = Dimensions[ren][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)
n = ren + I imn;

sinstep[t_] := 20 - (5/6 \[CapitalDelta]) Sin[\[Pi]/2 t]^2;
frames1 = Table[
  \[Phi]in = Table[E^(-((x)^2 + (y + sinstep[t])^2)/(2 (\[Lambda]0/5)^2)), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];(*Discretized source*)
  b = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
  M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k0^2]]; (*Operator on the left-hand side of the equation we want to solve*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  MatrixPlot[Transpose[(Re[(\[Phi]in + \[Phi]s)]/Max[Abs@Re[\[Phi]in + \[Phi]s][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]]])][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], ColorFunction -> ReMapC, DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, PlotLabel -> "Maxwell's fisheye: n=\!\(\*FractionBox[\(2\), \(1 + \
\*FractionBox[SuperscriptBox[\(r\), \(2\)], SuperscriptBox[\(R\), \(2\)]]\)]\)", LabelStyle -> {Black, Bold}, Epilog -> {White, Circle[{Round[dim/2 - (4 d)/\[Delta]], Round[dim/2 - (4 d)/\[Delta]]}, dim/3 - (1 d)/\[Delta]]}](*Plot everything*)
  , {t, 0, 1, 1/20}];

sinstep[t_] := 20 - (5/6 \[CapitalDelta]) Sin[\[Pi]/2 t]^2;
frames2 = Table[
  \[Phi]in = Table[E^(-((x - \[CapitalDelta]/3 Cos[\[Theta]])^2 + (y - \[CapitalDelta]/3 Sin[\[Theta]])^2)/(2 (\[Lambda]0/5)^2)), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];(*Discretized source*)
  b = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
  M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k0^2]]; (*Operator on the left-hand side of the equation we want to solve*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  MatrixPlot[Transpose[(Re[(\[Phi]in + \[Phi]s)]/Max[Abs@Re[\[Phi]in + \[Phi]s][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]]])][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], ColorFunction -> ReMapC, DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, PlotLabel -> "Maxwell's fisheye: n=\!\(\*FractionBox[\(2\), \(1 + \
\*FractionBox[SuperscriptBox[\(r\), \(2\)], SuperscriptBox[\(R\), \(2\)]]\)]\)", LabelStyle -> {Black, Bold}, Epilog -> {White, Circle[{Round[dim/2 - (4 d)/\[Delta]], Round[dim/2 - (4 d)/\[Delta]]}, dim/3 - (1 d)/\[Delta]]}](*Plot everything*)
  , {\[Theta], \[Pi]/2, 3/2 \[Pi], \[Pi]/20}];

sinstep[t_] := 1/3 \[CapitalDelta] + (\[CapitalDelta]/2 - \[CapitalDelta]/3) Sin[\[Pi]/2 t]^2;
frames3 = Table[
  \[Phi]in = Table[E^(-((x)^2 + (y + sinstep[t])^2)/(2 (\[Lambda]0/5)^2)), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];(*Discretized source*)
  b = -(Flatten[n]^2 - 1) k0^2 Flatten[\[Phi]in]; (*Right-hand side of the equation we want to solve*)
  M = L + DiagonalMatrix[SparseArray[Flatten[n]^2 k0^2]]; (*Operator on the left-hand side of the equation we want to solve*)
  \[Phi]s = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  MatrixPlot[ Transpose[(Re[(\[Phi]in + \[Phi]s)]/Max[Abs@Re[\[Phi]in + \[Phi]s][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]]])][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], ColorFunction -> ReMapC, DataReversed -> True, Frame -> False, PlotRange -> {-1, 1}, PlotLabel -> "Maxwell's fisheye: n=\!\(\*FractionBox[\(2\), \(1 + \
\*FractionBox[SuperscriptBox[\(r\), \(2\)], SuperscriptBox[\(R\), \(2\)]]\)]\)", LabelStyle -> {Black, Bold}, Epilog -> {White, Circle[{Round[dim/2 - (4 d)/\[Delta]], Round[dim/2 - (4 d)/\[Delta]]}, dim/3 - (1 d)/\[Delta]]}](*Plot everything*)
  , {t, 0, 1, 1/10}];

ListAnimate[Join[frames1, frames2, frames3]]