How to plot solution curves for an ODE with an implicit solution?

Hi, I have the given command

sol = Flatten[DSolve[y'[x] == y[x]^3/x^3 + y[x]/x + 1, y[x], x]]

which yields the solution

Solve[ArcTan[(-1 + (2 y[x])/x)/Sqrt[3]]/Sqrt[3] + 
   1/3 Log[1 + y[x]/x] - 1/6 Log[1 - y[x]/x + y[x]^2/x^2] == 
  C[1] + Log[x], y[x]]

When I try to plot solution curves for this, based on various C values, I do the following:

f[x_, y_] := sol
SolutionCurves[f_, xrange_, yrange_, init_, styl_ : {Red, Thick}] := 
StreamPlot[{1, f}, xrange, yrange, 
PlotRange -> {Take[xrange, -2], Take[yrange, -2]}, 
StreamPoints -> {init, Automatic, 10}, StreamScale -> None, 
StreamColorFunction -> False, StreamStyle -> styl]
SolutionCurves[
f[x, y], {x, -4, 4}, {y, -4, 
4}, {{1, -2}, {1.5, -1.01}, {2, 1}, {3, 2}}, {Blue, Medium}]
curveplot = %;

However, I get the error message

Solve::nsmet: This system cannot be solved with the methods available to Solve.
Part::partw: Part 1 of {} does not exist.

Any ideas on what I can do?

Thanks.

Cross-posted:

https://mathematica.stackexchange.com/questions/272907/error-message-on-plotting-a-solution-curve-for-an-ode-solved-with-dsolve