The stochastic differential equation,
$$ dX_t = (a_1(t)X_t + a_2(t))dt + b_2(t)dW_t, $$
has the general solution,
$$ X_t = \Phi_t\left( X_0 + \int_0^t a_2(s)\Phi_t^{-1}ds + \int_0^t b_2(s)\Phi_t^{-1}ds\right), $$
where
$$ \Phi_t = e^{\int_0^t a_1(s)ds}. $$
But where does this come from? How can we derive this? In a following post I'll describe how to find the more general solution when the SDE has multiplicative noise, but for now we're going to stick with the easier additive noise case.
To begin, let's go back to ODEs. If we had the equation,
$$ \frac{dX}{dt} - P(t)X = Q(t), $$
we can define the integrating factor
$$ I = e^{-\int_0^tP(s)ds}.$$
The reciprocal of the integrating factor solves the homogeneous equation, when the forcing term $Q(t)=0$. The point of this integrating factor is that it reformulates the problem into something easier. Notice that
$$\frac{d XI}{dt} = I\frac{dX}{dt} + X\frac{dI}{dt} = I\frac{dX}{dt} -PIX = I\left( \frac{dX}{dt} -PX\right) = IQ$$
So, integrating the far left and far right side of the above equation, we find that
$$XI = X(0) + \int_0^t I(s)Q(s)ds$$
since $I(0) = 1$. The solution the original ODE is
$$X = I^{-1}\left(X(0) + \int_0^t I(s)Q(s)ds\right)$$
Now we can apply this concept to the original SDE. This time the forcing terms are $a_2(t)dt$ and $b_2(t)dW_t$. So, the homogeneous equation is
$$dX^h_t = a_1(t)X^h_tdt$$
where the superscript $h$ denotes the homogeneous solution. Divide by $X_t^h$ and integrate to get
$$ln(X_t^h) = \int_0^t a_1(s)ds$$
Note that we can only do this because $X_t^h$ is a function of time only! See next post for more information on that point. Anyways, the homogeneous solution here is
$$X_t^h = e^{\int_0^t a_1(s)ds}$$
Let us call this homogeneous solution the fundamental solution and denote it by $X_t^h = \Phi_t$. It's reciprocal is the integrating factor we need. Playing the same trick as with ODEs, we evaluate $X_t\Phi_t^{-1}$ in the differential. But since we are dealing with stochastics, we need to be a little careful. The differential comes from taking a Taylor series expansion, so
$$X_{t+dt}(\Phi_{t+dt})^{-1} = X_t\Phi_t^{-1} + \frac{\partial X_t\Phi_t^{-1}}{\partial t}dt + \frac{\partial X_t\Phi_t^{-1}}{\partial W_t}dW_t + \frac{1}{2}\frac{\partial^2 X_t\Phi_t^{-1}}{\partial W_t^2}dt + h.o.t.$$
where the last term utilized Itô's Lemma that $dW_t^2 \to dt$. The term on the LHS and the first term on the RHS can be written simply as $dX_t\Phi_t^{-1}$. The fundamental solution is only a function of time, so we may rearrange the Taylor's series to get
$$dX_t\Phi_t^{-1} = X_t (-a_1\Phi_t^{-1})dt + \Phi_t^{-1}\left(\frac{\partial X_t}{\partial t} + \frac{1}{2}\frac{\partial^2 X_t}{\partial W_t^2}\right)dt + \Phi_t^{-1}\frac{\partial X_t}{\partial W_t}dW_t$$
which in turn reduces to
$$dX_t\Phi_t^{-1} = X_t (-a_1\Phi_t^{-1})dt + \Phi_t^{-1}dX_t.$$
Now, perhaps we could have simply written from the outset that $dX_t\Phi_t^{-1} = X_t d\Phi_t^{-1} + \Phi_t^{-1}dX_t$, and we would have been correct. But that would only be luck. We would need corrections when the SDE has multiplicative noise (next post).
Let's get back to the solution though. We have the fundamental solution $\Phi_t$ and we have that $dX_t\Phi_t^{-1}= X_t (-a_1(t)\Phi_t^{-1})dt + \Phi_t^{-1}dX_t$. Expanding $dX_t$ by using the original equation we were trying to solve reveals the coefficients to get
$$dX_t\Phi_t^{-1} = X_t (-a_1\Phi_t^{-1})dt + \Phi_t^{-1}\left((a_1X_t + a_2)dt + b_2dW_t\right)$$
Simplifying,
$$dX_t\Phi_t^{-1} = \Phi_t^{-1}\left(a_2dt + b_2dW_t\right)$$
Integrating,
$$X_t\Phi_t^{-1} = X_0 + \int_0^t (\Phi_s)^{-1}a_2(s)ds + \int_0^tb_2(s)(\Phi_s)^{-1}dW_s$$
since $\Phi_0^{-1} = 1$. Finally we have
$$X_t = \Phi_t\left(X_0 + \int_0^t (\Phi_s)^{-1}a_2(s)ds + \int_0^tb_2(s)(\Phi_s)^{-1}dW_s\right)$$
Let's end with an example. The Langevin equation is
$$dX_t = -\mu^2X_tdt + dW_t$$
Before we dive into the solution, it's always good practice to get a feel for what the equation is telling us. First, the system is not driven by some external function of time. Second, we have additive noise. This means that when $X_t$ is very large then the noise has little relative effect and when $X_t$ is very small the noise has much greater relative effect. Ignoring the noise, we have an exponential decay function. So we might expect that the function decays exponentially, but has noise (of the same magnitude throughout time) that jiggles the solution about the decaying function. Moreover, we are continually adding Gaussian noise, so we should expect (since Gaussian + Gaussian = Gaussian) that the stochastic process $X_t$ is itself Gaussian.
The solution is
$$X_t = X_0e^{-\mu^2t} + e^{-\mu^2t}\int_0^te^{\mu^2s}dW_s$$
Does this match our intuition? You tell me.
We can calculate mean and variance for this function rather easily. And since it's Gaussian, the mean and variance completely specify the density function for the process.