\[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/15; \[CapitalDelta] = 
 40*\[Lambda]0; (*Parameters for the grid*)
\[Sigma] = 
  10 \[Lambda]0; (*width of the gaussian beam*)

sourcef[x_, y_] := 
  E^(-(x^2/(2 \[Sigma]^2)))
    E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(I k0 y);
\[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*)

c1 = -\[CapitalDelta]/5; c2 = -\[CapitalDelta]/5;
ycenter = 
  Map[y0 /. # &, 
   FullSimplify[Solve[(x1)^2 + (y1 - y0)^2 == r^2, {y0}]][[All, 1, 
    All]]  ];
surface2[x_] := 
  Evaluate[Evaluate[((Sqrt[r^2 - (x)^2] + y0) /. {y0 -> 
          ycenter[[1]]}) /. {y1 -> -(\[CapitalDelta]/4), 
       x1 -> \[CapitalDelta]/2, 
       r -> (c1^2 + (c1 \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/
         16)/(2 (c1 + \[CapitalDelta]/4))} ] + 
    0.0 (Sin[3 x] + Sin[2 \[Pi] x])];
surface1[x_] := 
  Evaluate[((-Sqrt[r^2 - (x)^2] + y0 - 1) /. {y0 -> 
        ycenter[[2]]}) /. {y1 -> -(\[CapitalDelta]/4), 
     x1 -> \[CapitalDelta]/2, 
     r -> (c2^2 + (c2 \[CapitalDelta])/2 + (5 \[CapitalDelta]^2)/16)/(
      2 (c2 + \[CapitalDelta]/4))}];

ren = Table[
   If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 3, 
    1], {x, -\[CapitalDelta]/2, \[CapitalDelta]/
    2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/
    2, \[Delta]}];

stopstep[t_] := t (2 - t);
frames1 = Table[
   apodization = 
    Table[If[y < Re@Evaluate[surface2[x]] && y > Re@surface1[x], 
      1 (1 - E^(-(x^2/((200 - 190*stopstep[t]) \[Sigma]^2)))), 
      0], {x, -\[CapitalDelta]/2, \[CapitalDelta]/
      2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/
      2, \[Delta]}];
   n = ren + I imn + I apodization;
   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*)
   
   GraphicsRow[{
     ImageSubtract[
      ImageAdd[
       ArrayPlot[
        Transpose[(Abs[\[Phi]in + \[Phi]s]/
           Max[Abs[\[Phi]in + \[Phi]s]])^2][[(
          4 d)/\[Delta] ;; (-4 d)/\[Delta], (
          4 d)/\[Delta] ;; (-4 d)/\[Delta]]], 
        ColorFunction -> "AvocadoColors" , DataReversed -> True, 
        Frame -> False, PlotRange -> {0, 1}],
       ArrayPlot[Transpose@((ren - 1)/1) , DataReversed -> True , 
        ColorFunctionScaling -> False, ColorFunction -> GrayLevel, 
        Frame -> False]
       ],
      ArrayPlot[Transpose@(30*apodization) , DataReversed -> True , 
       ColorFunctionScaling -> False, ColorFunction -> GrayLevel, 
       Frame -> False]
      ]
     ,
     ListPlot[
      Transpose[(Abs[\[Phi]in + \[Phi]s]/
          Max[Abs[\[Phi]in + \[Phi]s]])^2][[(
         4 d)/\[Delta] ;; (-4 d)/\[Delta], (
         4 d)/\[Delta] ;; (-4 d)/\[Delta]]][[460, 
       271 - 100 ;; 271 + 100]], PlotRange -> All, Axes -> False, 
      Frame -> True, Joined -> True, FrameTicks -> None, 
      PlotStyle -> {Thick, Black}, ImageSize -> Large, 
      PlotLabel -> "Focus profile", 
      LabelStyle -> {Black, Bold, FontSize -> 16}]
     }]
   , {t, 0, 1, 1/100}];

ListAnimate[
 Join[frames1, Table[frames1[[-1]], {10}], Reverse@frames1, 
  Table[frames1[[1]], {5}] ] ]