Boundary value problem with a DiracDelta

@MichaelE2 gave the idea to use a shooting method, because NDSolve is only able to handle initial value problems involving DiracDelta

First solve the problem with a parametric slope f'[0]==fs0

F = ParametricNDSolveValue[{f''[x] + x f'[x] == x^2*DiracDelta[x - 1],
f[0] == 0, f'[0] == fs0}, f , {x, 0, 10}, fs0] 

Now choose fs0 to fullfill the second boundary condition f[10]==1

sol = FindRoot[F[fs0][10] == 1, {fs0, 1}]

Plot[F[fs0 /. sol][x], {x, 0, 10}, PlotRange -> All]

enter image description here


Replacing DiracDelta by its approximation in the weak topology helps:

s = NDSolve[{f''[x] + x f'[x] == x^2 *0.01/Pi/((x - 1)^2 + 0.01^2),
f[0] == 0, f[10] == 1}, f[x], {x, 0, 10}]
Plot[f[x] /. s, {x, 0, 10}, PlotRange -> All]

enter image description here


ClearAll[f, x];
psol = ParametricNDSolveValue[{f''[x] + x f'[x] == 
     x^2 DiracDelta[x - 1], f[0] == 0, f'[0] == p}, 
   f, {x, 0, 10}, {p}];

FindRoot[psol[p][10] == 1, {p, 0.2}]
bvpsol = psol[p] /. %;
(*  {p -> 0.274728}  *)

bvpsol[{0, 10}]
% - {0, 1}
(*
{0., 1.}
{0., -1.44329*10^-15}
*)