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).