Question: Method for numerical solving BVP of ODE with discontinuos right part - required

I try to solve numerically a boundary VP for ODE with different order of discontinuity of right part.

Say, the following BVP is given:

y''(x)+y'(x)+y(x)=F(x)

y(0)=1, y(2)=1

Let's use piecewise right part

F  := piecewise(x<=1, -x, x>1, 2*x+(x-1)^2)

The function

piecewise(x<=1, 1-x, x>1, (x-1)^2)

as obviuos, satisfies the BVP exclung the point x=1, where its 1st and 2nd derivatives are discontinuos.

Numerical solution

N0:=6:
As:=dsolve([diff(y(x), x$2)+diff(y(x), x)+y(x)=F,  y(0)=1, y(2)=1], y(x), type=numeric, output = Array([seq(2.0*k/N0, k=0..N0)]), 'maxmesh'=500, 'abserr'=1e-3):

provides the solution essentially different to exact one described above:

But if to use the right part

F := piecewise(x<=1, x^2+x+2, x>1, -x^2+x)

for which the function

piecewise(x<=1, 1-x+x^2, x>1, -1+3*x-x^2)

satisfies the BVP excluding x=1, where this function has discontinuity of 2nd derivative only, the corresponding numerical solution is very similar to this exact solution:

This reason of the difference between these two cases is clear. In the first case both 1st and 2nd derivatives are discontiuos, while in the second one -- 1st derivative is contiuos.

I wonder, if there are numerical methods, implemeted in Maple, for numerical solution of the first type BVP with non-smooth right part?