I want to compare different SDE solvers against the analytical solution for the Ornstein-Uhlenbeck process. Currently I am doing this by solving numerically first (e.g. Euler-Maruyama), saving the noise, and calculating the integral solution, e.g.

(1)Here is the code I am using:

```
using DifferentialEquations
using Plots
γ = 34.0
D = 69.0
μ= 5.3
u₀=7.0
f(u,p,t) = -γ*(u - μ)
g(u,p,t) = sqrt(D)
dt = 0.01
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,span)
sol = solve(prob,EM(),dt=dt,save_noise=true)
# Analytical solution
ex = exp.(-γ*sol.t)
inNoise = [0; diff(sol.W.u[1:end-1])]
anaSol = u₀*ex + μ*(1 .- ex) + sqrt(D)*ex.*cumsum(inNoise.*exp.(γ*sol.t))
plot(sol.t,anaSol,label = "Ana")
plot!(sol,label = "EM()")
```

This is the result I want, but it involves saving the noise and calculating Weiner increments for the integration. I know that `SDEProblem`

can be used to add a defined analytical solution to an SDE, but I’m not sure the solution (1) can be cast in this form. Is there a more elegant way to calculate the analytical solution?