nu = 
        0.3 ; L = 
        100 ; 
nel = 20;  a = 1/nel
  
pos = Table[-1 + (2  i - 1)/nel, {i, 1, nel}]

T1[c_, 
          x_, y_] = ArcTan[x - c - a, y];
T2[c_, x_, y_] = ArcTan[x - c + a, y];
RR1[c_, x_, y_] = (x - c - a)^2 + 
          y^2;
RR2[c_, x_, y_] = (x - c + a)^2 + 
          y^2;
sxx[c_, x_, y_] = 
        N[-1/
              Pi (T1[c, x, y] - T2[c, x, y] + 
                y (x - c - a)/RR1[c, x, y] - 
                y (x - c + a)/
                    RR2[c, x, y])];
syy[c_, x_, y_] = 
        N[-1/
              Pi (T1[c, x, y] - T2[c, x, y] - 
                y (x - c - a)/RR1[c, x, y] + 
                y (x - c + a)/
                    RR2[c, x, y])];
sxy[c_, x_, y_] = 
        N[-y^2/
              Pi (1/RR1[c, x, y]^2 - 
                1/RR2[c, x, y]^2)];
ux[c_, x_, y_] = 
        N[-1/(2  
                  Pi) ((1 - 
                      2  nu) ((x - c - a) 
                        T1[c, x, y] - (x - c + a) T2[c, x, y] - 
                      Pi a) + (1 - nu) 
                  y Log[RR1[c, x, y]/
                      RR2[c, x, y]]) ];
uy[c_, x_, y_] = 
        N[1/(2  Pi) (
-(1 - 2  nu) 
                  y (T1[c, x, y] - 
                      T2[c, x, y])
 + (1 - 
                      nu)
(  (x - c - a) 
                        Log[RR1[c, x, y]] - (x - c + a) 
                        Log[RR2[c, x, y]] + (L - c + a) 
                        Log[(L - c + a)^2] - (L - c - a) 
                        Log[(L - c - a)^2])  ) ];

(* 
      Fill the B matrix of influence coefficients *) 
(bmat = 
            Table[Table[ 
                uy[pos[[i]], pos[[j]], 0.0], {i, 1, nel}, {j, 1, 
                  nel}]]) // MatrixForm;
punch = Table[1, {i, 1, nel}];

pres = 
        LinearSolve[bmat, punch];

(* 
      Display the solution for contact pressure *) 
distrib = 
        Thread[{pos, pres}];
p1 = ListPlot[distrib, PlotStyle -> PointSize[0.02]]

(* 
      Introduce Ledge function for displaying the piecewise solution *) 
Ledge[c_, x_, d_] = 
        UnitStep[x - c + d] - 
          UnitStep[x - c - d];
Plot[Ledge[0, x, a], {x, -0.1, 0.1}]

(* 
      Plot the solution using the Ledge function *) 
Sol[x_] = 
        Sum[pres[[i]] Ledge[pos[[i]], x, a], {i, 1, 
            nel}];
p2 = Plot[Sol[x], {x, -1, 1}]

Show[{p1, p2}, 
      AxesOrigin -> {0. , 0. }]

(* 
      Evaluate displacement profile outside the punch *) 
points = 
        Table[1  + 10*(2  i - 1)/nel, {i, 1, 
            nel}];
(bbmat = 
            Table[Table[ 
                uy[points[[i]], pos[[j]], 0.0], {i, 1, 
                  Length[points]}, {j, 1, nel}]]) // 
        MatrixForm;
disp = Dot[bbmat, pres]
p3 = 
      ListPlot[Thread[{points, -disp}], 
        PlotStyle -> PointSize[0.02]]

rsyy[x_, y_] = 
        Dot[syy[pos, x, y], pres];
rsxx[x_, y_] = 
        Dot[sxx[pos, x, y], pres];
rsxy[x_, y_] = 
        Dot[sxy[pos, x, y], pres];
rux[x_, y_] = Dot[ux[pos, x, y], pres];
ruy[x_, y_] = Dot[uy[pos, x, y], pres];
dep = 0.01;
Plot[{rsyy[x, dep], rsxx[x, dep]}, {x, -2, 2}]

ContourPlot[rsyy[x, y], {x, -1.5, 1.5}, {y, 0.01, 3}, 
      PlotPoints -> 50]