How do I solve the first order nonlinear ODE? $$\dfrac{dy}{dx} = \dfrac{(1+ky-y/x)}{(1-y^2)}$$
$$y(0) = 0$$
I have tried using change of variables such as $x=y^2z$, or $x=yz$, however, I am still not able to separate the equation. Any thoughts on whether i should consider a different change of variable?
 A: This equation is not integrable for any real $k$.
Consider the planar system equivalent:
$$\begin{cases}x'&=&x(1-y^2)\\
y'&=&x+kxy-y\end{cases}\tag{1}$$
A solution to the original ODE is equivalent to a first integral for the system. By the theorem 5.2 in Goriely's "Integrability and nonintegrability of dynamical systems", this system has a first integral around any equilibrium point only if the linear eigenvalues of the Jacobian matrix are in "resonance", that is: for each equilibrium $p_j$, the eigenvalues $a_j$ and $b_j$ of the Jacobian matrix in $p_j$ must satisfy: $$ma_j+nb_j=0,$$
for $m$ and $n$ positive integers.
This system has equilibrium points: $p_1=(0,0)$, $p_2=(\frac{1}{k+1},1)$ and $p_3=(\frac{1}{k-1},-1)$. If you take the Jacobian in $p_2$, the eigenvalues are
$$a_2=\frac{-1+\sqrt{-8k^2-16k-7}}{2(k+1)}$$ and $$b_2=\frac{-1-\sqrt{-8k^2-16k-7}}{2(k+1)}.$$
So we need $m$ and $n$ that
$$m\frac{-1+\sqrt{-8k^2-16k-7}}{2(k+1)} +n\frac{-1-\sqrt{-8k^2-16k-7}}{2(k+1)}=0.$$
After some manipulations we get: $$8k^2+16k+7+(\frac{m+n}{m-n})^2=0,$$
wich has no real solutions.

Fig. 1: The vector field induced by (1) in the case $k=2$ displaying equilibrium points $(1/3,1)$ and $(1,-1)$ (figure included by JeanMarie).
Here is the Matlab program that has been used to create the plot:

axis([-2,2,-2,2]);hold on;
[x,y]=meshgrid(-2:0.1:2);
k=2;
u=x.*(1-y.^2);v=x+k*x.*y-y;
n=sqrt(u.^2+v.^2);
u=u./n;v=v./n;
quiver(x,y,u,v,1);


