(* This segment of code implements (3.1) in the paper *)
(* Note that we use 256 digits of accuracy for the Legendre functions.  This value can be adjusted. *)
g[n_]:=N[If[k2<4,((-1)^n/(2 π I))( LegendreQ[n-1/2,1-(4-k2)^2/8] - (1/2) π I LegendreP[n-1/2,1-(4-k2)^2/8]),
((-1)^n/(2 π I))( LegendreQ[n-1/2,1-(4-k2)^2/8] + (1/2) π I LegendreP[n-1/2,1-(4-k2)^2/8])],256]

(* k2 is the value of k squared, the same k squared on the left-hand side of (1.3a) *)
k2=1/2;

(* This next segment of code implements (3.2) in the paper *)
ClearAll[green,gg,ap,uap,u];
t = k2 - 4;
γ = -2;
green[0,0]=g[0];
For[m=1,m≤220,m+=1,
For[n=0,n≤m,n+=1,
If[(m==1) && (n==0),
green[1,0]=(2 t green[0,0] - 2)/(4 γ);green[0,1]=green[1,0],
If[m==n,
green[m,n]=g[n],
If[(m≠1) && (n==0),
green[m,0]=(2 t green[m-1,0] - γ green[m-2,0] - 2 γ green[m-1,1])/γ;green[0,m]=green[m,0],
If[m==(n+1),
green[m,n]=(2 t green[m-1,m-1] - 2 γ green[m-1,m-2])/(2 γ);
green[n,m]=green[m,n],
green[m,n]=(2 t green[m-1,n] - γ green[m-2,n] - γ green[m-1,n+1] - γ green[m-1,n-1])/γ;green[n,m]=green[m,n]
]
]
]
]
]
];
(* "gg" is the method of images Green's function *)
gg[i_,j_] := -green[Abs[i-1],Abs[j]]+green[Abs[i+1],Abs[j]];

(* "ap" specifies the size of the aperture *)
ap=10;

(* "f" specifies the function on the aperture, corresponding to (1.3b) in the paper *)
f[j_] := Piecewise[{
{1,Abs[j]≤ap},
{0,Abs[j]>ap}}];

(* "u" is the exact solution a la (2.22) in the paper *)
u[p_,q_] := If[p≥1,Sum[f[j]*gg[p,q-j],{j,-ap,ap}],uap[q]]

(* To 16 digits, this is the solution at i=100, i.e., 100 units to the right of the aperture, from j=-100 to j=100 *)
Table[N[Abs[u[100,j]],16],{j,-100,100}]

(* Here is a plot of the solution *)
ListPlot[Table[N[Abs[u[100,j]],16],{j,-100,100}]]

(* This function can be used to verify that the residual of the numerically computed solution is indeed zero *)test[i_,j_] := u[i,j+1] + u[i,j-1] + u[i+1,j] + u[i-1,j] - 4 u[i,j] + k2 u[i,j]

N[test[100,100],256]