diff --git a/Project.toml b/Project.toml
index e37edd243..254cf1ea5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FiniteHorizonGramians = "b59a298d-d283-4a37-9369-85a9f9a111a5"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 GaussianDistributions = "43dcc890-d446-5863-8d1a-14597580bb8d"
diff --git a/src/ProbNumDiffEq.jl b/src/ProbNumDiffEq.jl
index f8fc33edb..67c93476b 100644
--- a/src/ProbNumDiffEq.jl
+++ b/src/ProbNumDiffEq.jl
@@ -29,6 +29,7 @@ using Octavian
 using FastGaussQuadrature
 import Kronecker
 using ArrayAllocators
+using FiniteHorizonGramians
 
 @reexport using GaussianDistributions
 using GaussianDistributions: logpdf
diff --git a/src/priors/ltisde.jl b/src/priors/ltisde.jl
index 30ba4efda..a8f5e299a 100644
--- a/src/priors/ltisde.jl
+++ b/src/priors/ltisde.jl
@@ -37,7 +37,7 @@ function matrix_fraction_decomposition(
     return A, Q
 end
 
-function discretize_sqrt(sde::LTISDE, dt::Real)
+function discretize_sqrt_with_quadraturetrick(sde::LTISDE, dt::Real)
     F, L = drift(sde), dispersion(sde)
 
     D = size(F, 1)
@@ -71,3 +71,13 @@ function discretize_sqrt(sde::LTISDE, dt::Real)
 
     return Ah, Qh_R
 end
+
+
+function discretize_sqrt(sde::LTISDE, dt::Real)
+    Fdt, Ldt = dt*drift(sde), sqrt(dt)*dispersion(sde)
+
+    method = FiniteHorizonGramians.ExpAndGram{eltype(Fdt),13}()
+    Ah, Qh_R = FiniteHorizonGramians.exp_and_gram_chol(Fdt, Ldt, method)
+
+    return Ah, Qh_R
+end