\[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] = 50*\[Lambda]0; (*Parameters for the grid*) \[Sigma] = 10 \[Lambda]0; (*width of the gaussian beam*)
sourcef[x_, y_] :=E^(-((x + \[CapitalDelta]/8)^2 + (y + \[CapitalDelta]/3)^2)/(2 (\[Lambda]0/5)^2));
\[Phi]in = Table[Chop[sourcef[x, y]], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*Discretized source*)
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[\[Phi]in][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)
ReMapC[x_] := RGBColor[(2 x - 1) UnitStep[x - 0.5], 0, (1 - 2 x) UnitStep[0.5 - x]];

frames1 = Table[
  ren = Clip[
    Table[If[-\[Lambda]0 - \[CapitalDelta]/8 < x < \[Lambda]0 - \[CapitalDelta]/8, \[Alpha], 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {1, 2}];
  n = ren + I imn;
  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*)
ImageAdd[
   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}]
   ,
   ArrayPlot[Transpose[Re[(n - 1)/5]] [[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ]
  , {\[Alpha], 1, 2, 1/10}]

frames2 = 
 Table[ren = Clip[Table[If[-\[Lambda]0 - \[CapitalDelta]/8 < x < \[Lambda]0 - \[CapitalDelta]/8, 2, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}] + Table[If[\[CapitalDelta]/2 - 6*\[Lambda]0 < (x - \[CapitalDelta]/8 - \[Lambda]0/4 + \[CapitalDelta]/8 + (-(\[CapitalDelta]/1.7) (t - 1)^4))^2 + (y)^2 < \[CapitalDelta]/2 + 6*\[Lambda]0, 1, 0], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {1, 2}];
  n = ren + I imn;
  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*)
  ImageAdd[
   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}]
   ,
   ArrayPlot[Transpose[Re[(n - 1)/5]] [[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ]
  , {t, 0, 1, 1/10}]

frames3 = 
 Table[ren = Clip[Table[If[-\[Lambda]0 - \[CapitalDelta]/8 < x < \[Lambda]0 - \[CapitalDelta]/8, 2, 1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}] + Table[If[\[CapitalDelta]/2 - 6*\[Lambda]0 < (x - \[CapitalDelta]/8 - \[Lambda]0/4 + \[CapitalDelta]/8)^2 + (y)^2 < \[CapitalDelta]/2 + 6*\[Lambda]0, 1, 0], {x, -\[CapitalDelta]/2, \[CapitalDelta]/       2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}] + Table[If[-\[Lambda]0 + \[CapitalDelta]/4 + \[Lambda]0/2 - \[CapitalDelta]/8 + (\[CapitalDelta]/1.7 (t - 1)^4) < x < \[Lambda]0 + \[CapitalDelta]/4 + \[Lambda]0/2 - \[CapitalDelta]/8 + (\[CapitalDelta]/1.7 (t - 1)^4), 1, 0], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {1, 2}];
  n = ren + I imn;
  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*)
  ImageAdd[
   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}]
   ,
   ArrayPlot[Transpose[Re[(n - 1)/5]] [[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ]
  , {t, 0, 1, 1/10}]

ListAnimate[Join[
  Table[frames1[[1]], {5}], frames1,
  Table[frames2[[1]], {5}], frames2,
  Table[frames3[[1]], {5}], frames3, Table[frames3[[-1]], {20}]
  ], ImageSize -> Medium]