How to iterate all the points on curves

Question

I plot curves by

ContourPlot[x^2 y + x y^2 == 1, {x, -10, 10}, {y, -10, 10}]

It shows three curves on the (x,y) plane, but now my question is how to iterate all the points on the curves, the iterative equation is complex (nonlinear), for example, x = x+y*x, and y = x*y+y^2.

I have an idea, solve x^2 y + x y^2 == 1 and get (x,y), then plot it by ParametricPlot one by one, but it will be tedious (as it is nonlinear, so solve it is difficult sometimes). Is there any better way?

Answer 1

Update: An alternative approach using BSplineFunction:

ClearAll[tr, iter]

tr = {# + #2 #, # #2 + #2^2} & @@@ # &;

iter = Normal[#] /. Line -> (Line[tr[BSplineFunction[#] /@ Subdivide[50]]] &) &;

Examples:

cp = ContourPlot[x^2 y + x y^2 == 1, {x, -10, 10}, {y, -10, 10}];

NestList[iter, cp, 3]

NestList[Show[iter @ #, PlotRange -> All] &, cp, 3]

Original answer:

ClearAll[tr, iterate]

tr = {# + #2 #, # #2 + #2^2} & @@@ # &;
iterate = MapAt[tr, #, {1, 1, 1}] &;

Examples:

cp = ContourPlot[x^2 y + x y^2 == 1, {x, -10, 10}, {y, -10, 10}];

NestList[iterate, cp, 3]

NestList[Show[iterate @ #, PlotRange -> All] &, cp, 3]

Note: In version 11.3 replace {1, 1, 1} with {1, 1} in MapAt[...].

Answer 2

Clear["Global`*"]

There are two solutions

sol = Solve[x^2 y + x y^2 == 1, y] // Simplify

(* {{y -> -((x^(3/2) + Sqrt[4 + x^3])/(2 Sqrt[x]))}, {y -> -(x/2) + Sqrt[
    4 + x^3]/(2 Sqrt[x])}} *)

Plot[Evaluate[y /. sol], {x, -10, 10},
 Frame -> True,
 PlotLegends -> Placed["Expressions", {.75, .75}]]