\[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;(*size of a pixel*) \[CapitalDelta] = 30*\[Lambda]0; (*size of the grid*)
\[Sigma] = 3 \[Lambda]0; (*width of the source*)
d = \[Lambda]0/1; (*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[imn][[1]];
L = -1/\[Delta]^2*KirchhoffMatrix[GridGraph[{dim, dim}]]; (*Discretized Laplacian*)

ren0 = 1.; (*background refractive index*)
sourcef[x_, y_] := E^(-(x^2/(2 \[Sigma]^2))) E^(-((y + \[CapitalDelta]/2)^2/(2 (\[Lambda]0/2)^2))) E^(I k0 y); (*source founction*)
\[Phi]in = Table[Chop[sourcef[x, y] ], {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*input*)

(*a first run without nonlinearity to find a good normalization*)
ren = Table[
   ren0, {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}];
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]initial = Partition[LinearSolve[M, b], dim]; (*Solve the linear system*)
norm = Max[Abs[\[Phi]initial[[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]] ]^2]/2;

n2min = 0; (*minimum value of the nonlinear refractive index*)
n2max = N[1/35] (*maximum value of the nonlinear refractive index*)
\[CapitalDelta]n = Table[0, {x, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}, {y, -\[CapitalDelta]/2, \[CapitalDelta]/2, \[Delta]}]; (*initialize the change in refractive index to 0*)

n2 = n2min;
n2step = (n2max - n2min)/100;

evo = Reap[For[n2 = n2step, n2 <= n2max, n2 = n2 + n2step,
      \[Phi]old = \[Phi]initial; (*initialization*)
      \[Phi] = 2*\[Phi]initial; (*just to make the first iteration always run. It will be overwritten immediately*)
      While[
       Max[Abs[\[Phi] - \[Phi]old]] > 10^-3 (*stop iterating after the change per step get small*),
       \[Phi]old = \[Phi];
       intensity = (Abs[\[Phi]old]^2/norm)[[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]];
       \[CapitalDelta]n = n2*ArrayPad[intensity, (4 d)/\[Delta] - 1];
       n = ren + \[CapitalDelta]n + I imn; (*new refractive index profile*)
       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*)
       ];
      Sow[{\[CapitalDelta]n, \[Phi]}];
      ];][[2, 1]];

frames = Table[
   Grid[{{ArrayPlot[ evo[[j, 1]][[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], PlotRange -> {0, 0.05}, ClippingStyle -> Black, PlotLabel -> Style["Refractive index", Black, Bold]],       ArrayPlot[(Abs[evo[[j, 2]]]^2/norm)[[(4 d)/\[Delta] ;; (-4 d)/\[Delta], (4 d)/\[Delta] ;; (-4 d)/\[Delta]]], ColorFunction -> "AvocadoColors", PlotRange -> {0, 1.5}, PlotLabel -> Style["Intensity", Black, Bold]]}}]
   , {j, 1, Dimensions[evo][[1]] }];

ListAnimate[frames]