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: (I have added the corresponding analytic solution plot)

(*Boundary elements set up and material properties*) nb = nm = 20; nd=ne=nb+1; 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, 2/nb}];  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[-1, {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 = 6; 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]; plot1 = Plot[Sij . py, {x, 0, 3},PlotStyle -> Blue] AnalyticUy[h_] := -(1 - ((Log[h + Sqrt[(h^2) - 1]])/Log[2])) plot2 = Plot[AnalyticUy[h], {h, 1, 3}, PlotStyle -> {Dashed, Black}] plot3 = Plot[-1, {x, 0, 1}, PlotStyle -> {Dashed, Black}] Show[plot1, plot2, plot3] 

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