Solution (source code)

= Solution

If $X_0\sim N(0,(2\lambda)^{-1})$ independently of $B$, then
$$
\operatorname{Var}(X_t)
=e^{-2\lambda t}\frac1{2\lambda}
+\frac{1-e^{-2\lambda t}}{2\lambda}
=\frac1{2\lambda}.
$$
Thus $X_t\sim N(0,(2\lambda)^{-1})$ for every $t$. For $0<s<t$, the Markov decomposition
$$
X_t=e^{-\lambda(t-s)}X_s
+\int_s^te^{-\lambda(t-r)}\,dB_r
$$
has an increment independent of $X_s$, and hence
$$
\operatorname{Cov}(X_t,X_s)
=\frac{e^{-\lambda(t-s)}}{2\lambda}.
$$
This is the stationary Ornstein-Uhlenbeck covariance.