I can follow the rationale for integrating factors well enough, but they still feel like voodoo to me.
Every single description of integrating factors I've seen (and I've seen quite a few, including these) skips the question of what integrating factors are, and instead adopts a strictly heuristic approach: integrating factors are presented as a trick for solving certain differential equation problems, and the fact that "they work" is all one needs to know about them1.
Is there a way to present integrating factors in more fundamental terms, independent of their role in the solution of differential equations?
Given a general understanding of what a differential equation says, and before even considering the problem of "solving" such an equation, how would integrating factors enter into the picture?
Alternatively, what would be a plausible line of reasoning that would lead someone (who didn't already know about integrating factors) to consider them in the first place, either in the context of solving a differential equation or otherwise?
1 The appearance of integrating factors in formal mathematical discourse often resembles the hiring of some unsavory characters to take care of a nasty job. Sure, they are an embarrassment to the usual mathematical decorum, but they are put up with out of necessity: they get the job done. The embarrassment is made only worse by the fact that these "integrating factors" often show up in the company of so-called "inexact differentials", which are even more unsavory, if anything.
If you're comfortable with differential forms, here's the story as best I understand it.
The essential point is that a function $f$ of a single variable can be viewed as defining a path $\gamma(t) = (t,f(t))$ in the plane $\mathbb{R}^2$, so that if you can nicely rewrite a differential equation for $f$ in terms of $\gamma$, you can now use tools from differential geometry to get a better handle on your original equation. Indeed, consider the general first order linear differential equation for an unknown function $f$: $$ a(t)f^\prime(t) + b(t)f(t) + c(t) = 0. $$ If you define the path $\gamma$ in $\mathbb{R}^2$ by $\gamma(t) = (x(t),y(t)) := (t,f(t))$, then you can rewrite your differential equation as $$ \gamma^\ast\omega = 0 $$ for the $1$-form $$ \omega := (b(x)y+c(x))dx+a(x)dy, $$ since $$ \gamma^\ast\omega = (b(x(t))y(t)+c(x(t)))x^\prime(t)dt + a(x(t))y^\prime(t)dt = (a(t)f^\prime(t) + b(t)f(t) + c(t))dt. $$
Now, if $\omega$ is already exact, i.e., $\omega = dF$ for some scalar function $F$, then $$ 0 = \gamma^\ast \omega = \gamma^\ast dF = d(F \circ \gamma), $$ so that $F \circ \gamma$ is constant, and hence $F(t,f(t)) = C$ gives an implicit solution to the original equation. However, for $\omega$ to be exact, it must be closed, i.e., $d\omega = 0$, which in this concrete case yields $$ 0 = d\omega = (a^\prime(x)-b(x)) dx \wedge dy, $$ or equivalently, $$ a^\prime(x) = b(x). $$
So, what do you do if $\omega$ isn't even closed, i.e., if $a^\prime(x) \neq b(x)$? Well, you can try to find a nowhere-vanishing function $\mu$, an integrating factor, such that $\mu\omega$ is closed, i.e., $$ 0 = d(\mu \omega ) = d\mu \wedge \omega - \mu d\omega. $$ If you can find such a function $\mu$, then $\gamma^\ast \omega = 0$ if and only if $\gamma^\ast (\mu \omega) = 0$ where $\mu \omega$ is closed. If moreover, $\mu \omega$ is defined on a region in the plane over which every closed $1$-form is exact (e.g., your region is contractible), then $\mu \omega$ is exact, i.e., $\mu \omega = dG$ for some scalar function $G$, and $G(t,f(t))=C$ implicitly defines the solutions of your equation.
Now, in our concrete case, $$ 0 = d\mu \wedge \omega - \mu d\omega\\ = (\mu_x(x,y) dx + \mu_y(x,y) dy) \wedge ((b(x)y+c(x))dx+a(x)dy) - \mu(x,y)(a^\prime(x)-b(x))dx \wedge dy\\ = (a(x)\mu_x(x,y) - (b(x)y+c(x))\mu_y(x,y) - (a^\prime(x)-b(x))\mu(x,y))dx \wedge dy, $$ or equivalently, $$ a(x)\mu_x(x,y) - (b(x)y+c(x))\mu_y(x,y) - (a^\prime(x)-b(x))\mu(x,y) = 0. $$ In the special case where $a(x) = 1$, i.e., your original ODE was $f^\prime(t) + b(t)f(t) + c(t) = 0$, if you assume that $\mu = \mu(x)$, so that $\mu_y = 0$, then this reduces to $$ \mu^\prime(x) - b(x)\mu(x) = 0, $$ yielding the usual integrating factor $$ \mu(x) = \mu(x_0)e^{\int_{x_0}^x b(s)ds} $$ from the textbooks.
Finally, let me just point out that this story can be told in considerable generality. Let $M$ be a smooth manifold (e.g., $M$ an open subset of $\mathbb{R}^2$, as in the usual ODE textbooks), let $\omega \in \Omega^1(M)$ be a $1$-form, and consider the first order differential equation $$ \gamma^\ast \omega = 0 $$ for an unknown smooth path $\gamma : (a,b) \to M$. If $\omega$ is already exact, i.e., $\omega = dF$ for some $F \in C^\infty(M)$, then $$ 0 = \gamma^\ast \omega = \gamma^\ast dF =d(F \circ \gamma), $$ so that $F \circ \gamma$ is constant, and hence $F(\gamma(t)) = C$ gives an implicit equation for $\gamma$. For $\omega$ to be exact, it must be closed, i.e., $d\omega = 0$. If it isn't even closed, you can try to find a nowhere vanishing scalar function $\mu \in C^\infty(M)$, your integrating factor for $\omega$, such that $\mu\omega$ is closed, i.e., $$ 0 = d(\mu \omega) = d\mu \wedge \omega - \mu d\omega. $$ Then, if $H^1(M,\mathbb{R}) = 0$, so that every closed $1$-form is exact, you can find a function $G \in C^\infty(M)$ such that $\mu\omega = dG$, in which case $G(\gamma(t)) = C$ defines an implicit equation for $\gamma$. In particular, finding $\mu$ involves solving the potentially very pesky PDE $$ d\mu \wedge \omega - \mu d\omega = 0, $$ and all the "voodoo" consists of various special cases (traditionally for $M$ an open subset of $\mathbb{R}^2$) where this PDE is actually tractable.