Can NDSolve handle discontinuous data?

You can force NDSolve to use the finite element method:

uif = NDSolveValue[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
    u'[1] == 0}, u, {x, 0, 1}, 
   Method -> {"PDEDiscretization" -> "FiniteElement"}];

Possibly, this might be the default.

In any case if you compare to the analytical solution:

aif = DSolveValue[{Rationalize[-u''[x] == UnitBox[(x - 0.5) 0.5/0.2]],
    u[0] == 0, u'[1] == 0}, u, {x, 0, 1}];

The result looks good:

Plot[aif[x] - uif[x], {x, 0, 1}]

For a general case one should keep in mind, is that the quality of the solution will depend on the mesh having nodes or internal boundaries at the discontinuity. There is some documentation about this in the SolvingPDEwithFEM tutorial in the variable coefficients section and probably the mesh generation tutorial is also of interest.

Another option is to switch off "DiscontinuityProcessing"

NDSolve[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, 
  u'[1] == 0}, u, {x, 0, 1}, 
 Method -> {"DiscontinuityProcessing" -> False}]

to get the old behavior.

This is just an extended comment.

I'm not quite certain what's going on as we can easily implement a shooting method manually here.

(Shooting method, in short: parametrize in terms of u'[0] and very the parameter until u'[1] has the desired value.)

fun = ParametricNDSolveValue[{-u''[x] == UnitBox[(x - 0.5) 0.5/0.2], u[0] == 0, u'[0] == ud}, u, {x, 0, 1}, ud]

FindRoot[fun[ud]'[1], {ud, 1}]
(* {ud -> 0.4} *)

Plot[{UnitBox[(x - 0.5) 0.5/0.2], fun[0.4][x]}, {x, 0, 1}]

This gives the solution to you equation.

I've had difficulty before with UnitBox, although I've forgotten the exact context. UnitStep usually works better.

sol = NDSolve[{-u''[x] == UnitStep[x - 0.3] - UnitStep[x - 0.7], 
    u[0] == 0, u'[1] == 0}, u, {x, 0, 1}];

Plot[u[x] /. First[sol], {x, 0, 1}]