Return to Answer

The started as an extended comment but has morphed into an answer.

A numerical solution can be obtained directly by substituting r[u] = u^2 w[u], as suggested by the symbolic solution above. Then the ODE becomes

As already noted above, r[10^6] == -1 is not a valid proxy for r[Infinity] == -1. So, arbitrarily choose w[1] == 1. The r[u] is given by

This started as an extended comment but has morphed into an answer.

A numerical solution can be obtained directly by substituting r[u] = u^2 w[u], as suggested by the symbolic solution above. The ODE becomes

As already noted above, r[10^6] == -1 is not a valid proxy for r[Infinity] == -1. So, arbitrarily choose w[1] == 1. Then, r[u] is given by

The following is more an extended comment than an answer.

The started as an extended comment but has morphed into an answer.

Second Addendum

A numerical solution can be obtained directly by substituting r[u] = u^2 w[u], as suggested by the symbolic solution above. Then the ODE becomes

(Unevaluated[H[u] == u*D[r[u], u] - 2 r[u]] /. r[u] -> u^2 w[u]) // Simplify
(* H[u] == u^3 Derivative[1][w][u] *)

As already noted above, r[10^6] == -1 is not a valid proxy for r[Infinity] == -1. So, arbitrarily choose w[1] == 1. The r[u] is given by

sol = u^2 w[u] /. Flatten@NDSolve[{w'[u] == H[u]/u^3, w[1] == 1}, w, {u, 10^-6, 2}, 
    AccuracyGoal -> Infinity]

plus C[1] u^2 as obtained symbolically earlier.

Plot[sol, {u, 10^-6, 2}, ImageSize -> Large, AxesLabel -> {u, r}, 
    LabelStyle -> Directive[12, Bold, Black]]

Experimentation suggests that the double integral defining H[u] behaves a bit better at small u with options, Exclusions -> {0, 0}, AccuracyGoal -> Infinity, PrecisionGoal -> 6. Then, H[u]/u^3 scales roughly as u^(-6/7}, as can be seen from

The integral can, then, be performed (as requested by the OP in a comment below), although its accuracy is a bit uncertain.

NIntegrate[H[u]/u^3, {u, 0, 1}]
(* 21636.9 *)

Experimentation suggests that the double integral defining H[u] behaves a bit better at small u with options, Exclusions -> {0, 0}, AccuracyGoal -> Infinity, PrecisionGoal -> 6, MaxRecursion -> 12. Then, H[u]/u^3 scales roughly as u^(-6/7}, as can be seen from

The integral can, then, be performed without error messages (as requested by the OP in a comment below).

NIntegrate[H[u]/u^3, {u, 0, 1}, MaxRecursion -> 12]
(* 21636.9 *)