With modified mesh Needs["NDSolve`FEM`"]; additionalPoints = {{Ls/2 }} beamDomain =ToElementMesh[ImplicitRegion[0 <= x <= Ls, x], "IncludePoints" -> additionalPoints, "MeshOrder" -> 1] and load q[x_] := Piecewise[{{qs, 0 <= x <= lss}, {qss, lss < x <= Ls}}] we can confirm result derived by @BillWatts using FiniteElement, if we include correct boundary conditions `moment[Ls/2]==-(219375/32)`. solution = NDSolveValue[{{D[moment[x], {x, 2}] + q[x] == 0, moment[x] == -Es Is*D[w[x], {x, 2}]}, {w[0] == 0, DirichletCondition[w[x] == 0, x == Ls/2], w[Ls] == 0, moment[0] == 0, DirichletCondition[moment[x] == -(219375/32), x == Ls/2], moment[Ls] == 0}}, {w, moment}, Element[x, beamDomain] , Method -> {"FiniteElement"} ] GraphicsRow[{Plot[solution[[1]][x], {x, 0, Ls}], Plot[solution[[2]][x], {x, 0, Ls}]}] [![enter image description here][1]][1] [1]: https://i.sstatic.net/Kny1eUUG.png