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}]}]
The question remains how to get these inner boundary conditions directly without friendly support by classical beam theory (thanks @Bill Watts) ?
I tried, still without success, to find the NeumannValue of this problem.
Any ideas from FEM experts?