I am trying to numerically solve an integral equation in Starfield and Crouch textbook (Boundary Element Methods in Solid Mechanics: With Applications in Rock Mechanics and Geological Engineering), equation 3.2.5 (see code). After computing the density py I attempted to plot uy vs x. I observed some numerical instabilities in the region 0<=x<=1. I Increased my number of integration points,np and also reduced my boundary element size,xe, I still observe the same irregularity.I’m not sure why this is happening. Is there something I’m not doing right? I would appreciate any guidance. Here is my code:

(*Boundary elements set up and material properties*)
nb = 10; nd = 11; ne = 11; nm = 10; G = 1; v = 0.1; L = 1.25;a = 0.1; 
(*Input the coordinates of the of ends of boundary elements (xe,ye)*)
xe = Table[i, {i, -1, 1, 0.2}]; 
ye = Table[0, {i, 1, ne}]; 
(*Input the coordinates of the midpoints of boundary elements (xm,ym)*)

xm = ym = Table[0, {i, 1, nm}]; 
jb = If[j < ne, j + 1]; 
Do[xm[[j]] = (xe[[j]] + xe[[jb]])/2; 
  ym[[j]] = (ye[[j]] + ye[[jb]])/2, {j, 1, nm}];
bv = Table[-0.001, {i, 1, nm}]; 

(*Compute elements of Influence coefficients Bij and Sij*)

Sij = Bij = Table[0, {i, 1, nb}, {j, 1, nb}]; 
uy = (1/(2 G Pi)) (-2 (1 - v) (Log[Sqrt[(x - xi)^2 + y^2]] - 
    Log[L - xi]) + y^2/((x - xi)^2 + y^2))(*Equation 3.2.5*); 
Get["NumericalDifferentialEquationAnalysis`"]; 
np = 10; points = weights = Table[Null, {np}]; 
Do[points[[i]] = GaussianQuadratureWeights[np, -1, 1][[i, 1]], {i, 1, np}]
Do[weights[[i]] = GaussianQuadratureWeights[np, -1, 1][[i, 2]], {i, 1, np}]
GuassInt[f_, z_] := Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]
Do[xb = (1/2)*(xe[[jb]]*(1 - z) + xe[[j]]*(1 + z)); yb = (1/2)*(ye[[jb]]*(1 -z) + ye[[j]]*(1 + z)); 
Do[Bij[[i, j]] = GuassInt[uy /. {x -> xm[[i]], xi -> xb, y -> yb}, z]; Sij[[i, j]] = GuassInt[uy /. {x -> x, y -> yb, xi -> xb}, z], {i, 1, nb}], {j, 1, nb}]
py = LinearSolve[Bij, bv];
Plot[Sij . py/0.001, {x, 0, 3}]

Here is my plot (see the red ellipse region).