<< NDSolve`FEM`

λ = 0.53; k0 = 2 π/λ; R = λ;

mesh = ToElementMesh[FullRegion[2], {{0, R}, {0, R}},
   "MaxCellMeasure" -> 0.0005];

mesh["Wireframe"]

op =
  Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] -
   k0^2 {u[x, y], v[x, y], 0}]

pde = op == {0, 0};

Subscript[Γ, D] =
  DirichletCondition[{u[x, y] == 0., v[x, y] == Exp[I k0 x]}, True];

{us, vs} =
  NDSolveValue[{pde, Subscript[Γ, D]}, {u, v}, {x, y} ∈ mesh]

DensityPlot[Re[vs[x, y]], {x, y} ∈ mesh,
  ColorFunction -> "Rainbow",
  PlotLegends -> Automatic,
  PlotPoints -> 50,
  PlotRange -> All]

λ = 53/100; k0 = 2 π/λ; R = λ;

op = Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] -
    k0^2 {u[x, y], v[x, y], 0}];

pde = op == {0, 0};
bc = Function[{x, y}, {u[x, y] == 0, v[x, y] == Exp[I k0 x]}] @@@ {{0, y}, {R, y}, {x,
      0}, {x, R}} // Flatten;

domain = {0, R};
points = 50;
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoae isn't included in this code piece,
  please find it in the link above. *)
ptoa = pdetoae[{u, v}[x, y], {grid, grid}, difforder];
del = Most@Rest@# &;

ae = del /@ del@# & /@ ptoa@pde;
aebc = MapAt[del, ptoa@bc, List /@ Range@4];
{b, mat} = CoefficientArrays[{ae, aebc} // Flatten,
   Outer[#[#2, #3] &, {u, v}, grid, grid] // Flatten];
sollst = LinearSolve[-mat, N@b];
solmat = ArrayReshape[sollst, {2, points, points}];
{solu, solv} = ListInterpolation[#, {grid, grid}] & /@ solmat;

Outer[Plot3D[#@First[#2][x, y], {x, 0, R}, {y, 0, R}, PlotLabel -> #@Last@#2] &, {Re,
   Im}, {{solu, "u"}, {solv, "v"}}, 1] // Grid