p[x_] = 7.95775*10^-6/(6.25*10^-10 + x^2) /. x -> 10 - x

eq1 = D[M[x], {x, 2}] == p[x]

soln1 = NDSolve[{eq1, M[0] == 0, M[20] == 0}, M, {x, 0, 20}]

Plot[-M[x] /. soln1[[1]], {x, 0, 20}, PlotRange -> {0, 6}, FrameLabel -> {"x", "M[x]"}, PlotLabel -> "Bending moment", GridLines -> Automatic, Frame -> True]

V[x_] = -D[M[x] /. soln1[[1]], {x, 1}]

Plot[V[x], {x, 0, 20}, PlotRange -> {-0.75, 0.75}, FrameLabel -> {"x", "V[x]"}, PlotLabel -> "Shear force", GridLines -> Automatic, Frame -> True]