\begin {align*} y^{\prime } & =\frac {y-xf\left ( x^{2}+ay^{2}\right ) }{x+ayf\left ( x^{2}+ay^{2}\right ) }\\ & =\omega \left ( x,y\right ) \end {align*}
Using anstaz’s it is found that \begin {align*} \xi & =-ay\\ \eta & =x \end {align*}
Hence\begin {align} \frac {dx}{\xi } & =\frac {dy}{\eta }=dS\nonumber \\ \frac {dx}{-ay} & =\frac {dy}{x}=dS \tag {1} \end {align}
The first two give\[ \frac {dy}{dx}=\frac {x}{-ay}\] This is separable. Solving gives (taking one root)\[ y=\frac {\sqrt {a\left ( ac_{1}-x^{2}\right ) }}{a}\] Solving for \(c_{1}\) gives\[ c_{1}=\frac {x^{2}+ay^{2}}{a}\] Hence\[ R=\frac {x^{2}+ay^{2}}{a}\] \(S\) is found from either \(\frac {dy}{\eta }=dS\,\) or \(\frac {dx}{\xi }=dS\). Using \(\frac {dx}{-ay}=dS\) then\[ \frac {dx}{-ay}=dS \] But \(y=\frac {\sqrt {a\left ( ac_{1}-x^{2}\right ) }}{a}\). Hence\begin {align*} \frac {dx}{-a\frac {\sqrt {a\left ( ac_{1}-x^{2}\right ) }}{a}} & =dS\\ \frac {dx}{-\sqrt {a\left ( ac_{1}-x^{2}\right ) }} & =dS\\ -\frac {1}{\sqrt {a}}\arctan \left ( \frac {\sqrt {a}x}{\sqrt {c_{1}a^{2}-x^{2}a}}\right ) & =S\\ -\frac {1}{\sqrt {a}}\arctan \left ( \frac {\sqrt {a}x}{ay}\right ) & =S \end {align*}
Where constant of integration is set to zero. What is left is to find \(\frac {dS}{dR}\). This is given by\begin {equation} \frac {dS}{dR}=\frac {S_{x}+S_{y}\omega }{R_{x}+R_{y}\omega } \tag {2} \end {equation} But\begin {align*} R_{x} & =\frac {2x}{a}\\ R_{y} & =2y\\ S_{x} & =-\frac {y}{x^{2}y^{2}+a}\\ S_{y} & =-\frac {x}{a\left ( 1+\frac {x^{2}y^{2}}{a}\right ) } \end {align*}
Hence (2) becomes\[ \frac {dS}{dR}=\frac {-\frac {y}{x^{2}y^{2}+a}+\left ( -\frac {x}{a\left ( 1+\frac {x^{2}y^{2}}{a}\right ) }\right ) \omega }{\frac {2x}{a}+2y\omega }\] But \(R=\frac {x^{2}+ay^{2}}{a}\). The above becomes\[ \frac {dS}{dR}=\frac {-\frac {y}{aR}+\left ( -\frac {x}{a\left ( 1+\frac {x^{2}y^{2}}{a}\right ) }\right ) \omega }{\frac {2x}{a}+2y\omega }\] To finish. Another hard part of this Lie method is to convert back \(\,\frac {dS}{dR}=\frac {S_{x}+S_{y}\omega }{R_{x}+R_{y}\omega }\) so that the RHS is only a function of \(R\). Need to find a robust way to do this. This is now a weak point in my program as I have few ode’s that it can’t do it