Heyho,
I am using the separation-of-variables method for quite a while now, but what was always bothering me a bit, is why is it possible to do those operations. I'll give a concrete example (source Wikipedia):
$$\frac{dy}{dx}=xy^2 + x \Rightarrow \frac{dy}{1+y^2} = x \:dx \Rightarrow \int{\frac{dy}{1+y^2}} = \int{x \:dx} \Rightarrow \cdots$$
and so on. My problem lies in step 2. Why can I just treat the differential operator like a variable?
And besides all those "technical" answers: For a differential equation, once you have an alleged solution, you can easily check it. So even if you used strange, non-rigorous, mystical methods, that last step shows whether or not your solution is correct.
Maybe the best way to see what's going on is to write it as an "almost separation of variables", where you separate the variables except that you keep the derivative intact: $$\frac{\mathrm dy}{\mathrm dx}=xy^2 + x \implies \frac{1}{1+y^2}\frac{\mathrm dy}{\mathrm dx} = x$$ This is clearly a valid transformation. And it is also valid to integrate both sides over the same variable: $$\frac{1}{1+y^2}\frac{\mathrm dy}{\mathrm dx} = x \implies \int \frac{1}{1+y^2}\frac{\mathrm dy}{\mathrm dx}\,\mathrm dx = \int x\,\mathrm dx$$ Now on the integral of the left hand side, you can do a variable substitution from $x$ to $y$. This gives $$\int \frac{1}{1+y^2}\frac{\mathrm dy}{\mathrm dx}\,\mathrm dx = \int \frac{1}{1+y^2}\,\mathrm dy$$ and therefore you get $$\int \frac{1}{1+y^2}\,\mathrm dy = \int x\,\mathrm dx$$ But that's exactly what you get by doing separation of variables the usual way.
To split the derivative $\frac{dy}{dx}$ in $dy$ and $dx$ is just a convenient formal way to solve the equation. If you do not want to go far into differentials etc, you can justify this splitting by the chain rule. For example, your equation can be rewritten as $$ y'(x)=xy(x)^2 + x \Rightarrow \frac{1}{1+y(x)^2}\cdot y'(x) =x \Rightarrow \arctan(y(x))' = x \Rightarrow \dots $$ since $\arctan(y)=\int\frac{dy}{1+y^2}$. Compare to see that the last expression looks exactly like in your post.
The "recipe" you are unhappy with can be validated as follows:
After step 1 your ODE has the form $$f(y)\>dy=g(x)\>dx\ .$$ This is saying that the quantities $x$ and $y$, when considered as functions of a "hidden" variable $t$, satisfy $$f\bigl(y(t)\bigr)\>\dot y(t)\ \equiv\ g\bigl(x(t)\bigr)\>\dot x(t)\tag{1}$$ for all $t$. Integrating $(1)$ with respect to $t$ from $t=0$ to some chosen upper limit $T$ we obtain $$\int_0^T f\bigl(y(t)\bigr)\>\dot y(t)\>dt=\int_0^T g\bigl(x(t)\bigr)\>\dot x(t)\>dt\ .$$ Assume that $f$ and $g$ have known primitives $F$ and $G$. Then a simple substitution in $(2)$ shows that we have $$F\bigl(y(T)\bigr)-F\bigl(y(0)\bigr)=G\bigl(x(T)\bigr)-G\bigl(x(0)\bigr)$$ for all $T$ in a suitable neighborhood of $0$. But this is saying that for all such $T$ the values $y:=y(T)$ and $x:=x(T)$ are related by $$F(y)-F(y_0)=G(x)-G(x_0)\ .$$
$$\frac{dy}{dx}=x(y^2 + 1) \Rightarrow \frac{dy}{1+y^2} = x \:dx$$
One can write $$ f'(x) = x((f(x))^2+1) $$ and then $$ \frac{f'(x)}{1+f(x)^2} = x. $$ Hence $$ \int \frac 1 {1+f(x)^2}\Big(f'(x)\,dx\Big) = \int x\,dx $$ $$ \int \frac 1 {1+u^2}\, du = \int x\,dx. $$
In other words, you're using the chain rule.
One can also write $$ \frac{\Delta y}{\Delta x} \approx x(y^2+1) \tag 1 $$ $$ \frac{\Delta y}{1+y^2} \approx x\,\Delta x. \tag 2 $$ The limit of $(1)$ as $\Delta x\to0$ is the derivative $dy/dx$ and the limit as $\Delta x\to0$ of a sum of values of $(2)$ is the integral.
If we have $$\frac{dy(x)}{dx}=f(x)g(y(x)),$$ then the careful way of solving this equation (assuming $g(y) \neq 0$) is to divide by $g(y)$ in order to get
$$\frac{1}{g(y(x))} \frac{dy(x)}{dx}=f(x).$$
Replace $x$ with $t$ and integrate each side from an arbitrary lower bound $a$ to $x$:
$$\int_a^{x} \frac{1}{g(y(t))} \frac{dy(t)}{dt} dt= \int_a^{x} f(t) dt$$
Now note that on the LHS we can do the substitution $w=y(t)$, which gives $dw=y'(t) dt$ so that the above becomes
$$\int^{y(x)}_{y(a)} \frac{1}{g(w)} dw =\int_a^x f(t) dt.$$ Now be sloppy with notation and replace $t$ with $x$, $w$ with $y$ and write $y(x)$ as $x$ and $b=y(a)$ so that we have
$$\int^{y}_{b} \frac{1}{g(y)} dy =\int_a^x f(x) dx$$ or $$\int \frac{1}{g(y)} dy =\int f(x) dx + C.$$
You'll get the exact same thing as if you split up the $dy/dx$, even though it isn't rigorous, so why spend all that time doing it the long way?
Would you have a problem with this if you replaced "d"s with "delta"s (finite differences) and integrals with finite sums? If not, just take the limit for smaller "delta"s.
I personally see dx or dy as just finite small numbers. Doing a little Matlab numerical differentitation would probably make you more comfortable with this.
