\[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] = 30*\[Lambda]0; (*Parameters for the grid*)

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[imn][[1]];

ren = Clip[Total[Table[ RotateRight[DiskMatrix[8, dim], {RandomInteger[{-Round[dim/2], Round[dim/2] }], RandomInteger[{ -Round[dim/2]/4, 0}]}], {70}]], {0, 1}] + 1;
d = \[Lambda]0/2; (*typical scale of the absorbing layer*)
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)
n = ren + I imn;
optimizedpoint = Round[0.8 dim];

sourcelist = Table[x0, {x0, -\[CapitalDelta]/2 + 2, \[CapitalDelta]/2 - 2, 1}];
sourcedim = Dimensions[sourcelist][[1]];

phases = Table[
  \[Phi]in = Table[E^(-((x - sourcelist[[j]])^2/(2  (\[Lambda]0/2)^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^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] = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  
  Arg[\[Phi][[Round[dim/2] , optimizedpoint ]] ]
  , {j, 1, sourcedim, 1}];

frames = Table[
  \[Phi]in = Total@Table[ Table[E^(-((x - sourcelist[[j]])^2/(2  (\[Lambda]0/2)^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(-I*t* phases[[j]]), {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}], {j, 1, sourcedim}];(*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] = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
  ImageAdd[
   ArrayPlot[Transpose[(Abs[(\[Phi])]^2/Max[(Abs[\[Phi]]^2)[[(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 -> "AvocadoColors", DataReversed -> True, Frame -> False, PlotRange -> {0, 1}, LabelStyle -> {Black, Bold}, ColorFunctionScaling -> False, Epilog -> {Red, Thick, Table[ Line[{{(k - 1) (dim - (8 d)/\[Delta])/sourcedim, t*5*phases[[k]] + 10}, {k (dim - (8 d)/\[Delta])/sourcedim, t*5*phases[[k]] + 10}}], {k, 1, sourcedim}] }],
   ArrayPlot[Transpose@Re[(tmpn - 1)/5] , DataReversed -> True , ColorFunctionScaling -> False, ColorFunction -> GrayLevel, Frame -> False]
   ]
  , {t, 0, 1, 1/40}]

ListAnimate[ Join[Table[frames[[1]], 15], frames, Table[frames[[-1]], 15], Reverse[frames]] ]