eq1 = 1/5 (y[x] - 2 x y'[x]) == D[(x y'[x])/y[x] + x y[x]^3 D[D[x y'[x], x]/x, x], x]/x;
bc1 = y[x] - 2 x y'[x] == 0;
bc2 = y[x] + 4 x^2 y''[x] == 0;

solSeries = seriesDSolve[eq1, y, {x, 0, 5},
                         {y[0] -> a, y'[0] -> 0, y''[0] -> b, y'''[0] -> 0}]

newbclist = Thread[(Derivative[#][y][x0] ==
        (D[Normal@solSeries, {x, #}] /. x -> x0) & ) /@ Range[0, 3]]

x0 = 1/10^4; xMax = 10;
sol = ParametricNDSolveValue[{eq1, newbclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> sol[a, b] // Evaluate,
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]

bclist = Thread[(Derivative[#][y][x0] & /@ Range[0, 3]) == {a, 0, b, 0}]

solfake = ParametricNDSolveValue[{eq1, bclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> solfake[a, b] // Evaluate,
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]