Bifurcation of integral curves Consider the following first order ODE:
$$\frac{\operatorname{d}\!y}{\operatorname{d}\!x} = x^2 - y^2$$
Despite the fact that this ODE has a very simple expression, it is not solvable in terms of elementary functions. (We need the so-called Bessel function $J_u(z)$, where $u \in \mathbb{R}$ and $z \in \mathbb{C}$.)
I've used AutoGraph to plot the direction field and to plot several integral curves.

There is a clear separation of the plane. Your eye will no-doubt pick out two very heavily coloured curve-like regions. The underlying curves form a bifurcation set: choosing points on either side give qualitatively different integral curves through those points. 
Is there a general way to find an equation or a parametrisation for the bifurcation set? Or do we have to be able to solve the ODE explicitly?
If it isn't possible to find the bifurcation set explicitly, then is there any way to find other information, e.g. how many regions the bifurcation set separates the plane into? (In my example, the plane is separated into three regions.)
Addendum: Here's a plot to show that the lines $y = \pm x$ ($x^2-y^2=0$) have no local significance to the integral curves. They lines $y=\pm x$ do seem to be asymptotes for some of the integral curves.

 A: Solving of $\frac{dy}{dx}=x^2-y^2$
Let $y=\frac{F'}{F}$
$y'=\frac{F''}{F}-\frac{F'^2}{F^2}=x^2-\frac{F'^2}{F^2}$
$$F''-x^2 F=0$$
The general solution of this parabolic cylinder ODE is :
$$F(x)=c_1D_{-1/2}(\sqrt{2}x)+c_2D_{-1/2}(i\sqrt{2}x)$$
where $D_\nu(x)$ is the parabolic cylinder function. Alternatively, it can be expressed with Bessel functions of order $\pm\frac{1}{4}$
$$y(x)=\frac{F'}{F}=\frac{c_1\left( x D_{-1/2}(\sqrt{2}x)-\sqrt{2} D_{1/2}(\sqrt{2}x)\right) + c_2\left( -x D_{-1/2}(i\sqrt{2}x)-i\sqrt{2} D_{1/2}(i\sqrt{2}x)  \right) }{c_1 D_{-1/2}(\sqrt{2}x)+c_2D_{-1/2}(i\sqrt{2}x)}$$
Note that there is only one independant constant $C$ by simplification with ether $C=\frac{c_1}{c_2}$ or $C=\frac{c_2}{c_1}$
I fact, the purpose of is the asymptotic behaviour of $y(x)$.
The asmptotic expansion of the parabolic cylinder function is :
$$D_{-1/2}(X)=e^{-X^2/4}\sqrt{X}\left(1-\frac{3}{8}\frac{1}{X^2}+O\left(\frac{1}{X^3} \right)   \right)$$
Applying this formulas to the above equations $F(x)$ and $F'(x)$, then $y(x)=\frac{F'}{F}$ is an arduous task which leads to 
$$y=\pm \left(x-\frac{1}{2x}-\frac{3}{8x^2}+O\left(\frac{1}{x^3} \right)\right) $$
The asymptotes are $y=x$ and $y=-x$
Of course, one must no confuse the asymptotes with the asymptotic or "boudary" curves (figure below) : $y=x$ and $y=-x$ are asymptotes of both boundary curves and integral curves.

NOTE :
The parabolic cylinder functions are related to some Bessel functions. So, it is possible to express $y(x)$ on equivalent forms with ether one or the other of those two kind of special functions. The choice of the parabolic cylinder leads to simpler calculus as far as we are studying the asympotic behaviour because the asymptotic series expansion is easier. On the other hand, it involves complex coefficients $c_1$ , $c_2$. So, the formula with the parabolic cylinder functions is not pleasant to use if we want to compute particular values in the real range of $y(x)$. For this different purpose, the formula with Bessel functions is more convenient :
$$y(x)=\frac{\left( C_1\left( x^2 I_{-3/4}(x^2/2)+x^2 I_{5/4}(x^2/2)+I_{1/4}(x^2/2)\right) + C_2\left(  x^2 I_{3/4}(x^2/2)+x^2 I_{-5/4}(x^2/2)+I_{-1/4}(x^2/2)\right) \right) }{2 x \left(C_1 I_{1/4}(x^2/2)+C_2 I_{-1/4}(x^2/2) \right)}$$
The symbol $I$ denotes the modified Bessel functions of the first kind.
For example, the figure below shows the particular curves corresponding to $C_1=0$ and $C_2=0$ respectively. 

