A stochastic process $X(t)$ by definition is gaussian iff all its finite-dimensional joint probability density functions are multivariate gaussian. Namely iff given the times $(t_1,t_2,...,t_n)$ , the random variables $(X(t_1),X(t_2),...,X(t_n))$ are jointly gaussian (n is arbitrary).
My intent is understanding why the solution of the following simple stochastic differential equation is a gaussian process:
$$dX(t)=a(t)dt+b(t)dW(t)$$
where W(t) is a Wiener process with diffusion coefficient D and the initial condition is $X(t_0)=X_0$, where $X_0$ is a random variable that is independent of W(t). a(t) and b(t) are deterministic "regular" functions.
The explicit solution is:
$$X(t)=X_0+ \int_{t_0}^{t}a(s)ds+\int_{t_0}^{t}b(s)dW(s)$$
If the integral $\int_{t_0}^{t}b(t)dW(t)$ were a gaussian process (I'm not sure about this) and the initial condition were a gaussian random variable, than at fixed t, the random variable $X(t)$ would be the sum of 3 independent gaussian random variables ( $\int_{t_0}^{t}a(t)dt$ is in fact degenerate) and would be therefore gaussian. But what about the multiple-times joint densities of $(X(t_1),X(t_2),...,X(t_n))$?
You can get rid of the drift applying the Girsanov's theorem (see for instance Oksendal's introductory book on SDEs, or "Introduction to stochastic integration" by H.H. Kuo), so you end up with $X(t)$ being a Wiener integral (plus the initial condition, but since it's independent and Gaussian I will ignore it hereafter ).
As you may know the Wiener integrals are gaussian random variables, and by the linearity of the integral any linear combination of Wiener integrals will be in turn a Wiener integral and hence gaussian. This allows to conclude that variables $X(t_1),...X(t_n)$ are jointly gaussian.