diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl
index d7c75bb..e7ce61e 100644
--- a/src/SpinShuttling.jl
+++ b/src/SpinShuttling.jl
@@ -50,7 +50,7 @@ function Base.show(io::IO, model::ShuttlingModel)
     println(io, "Process Time: T=$(model.T)")
     println(io, "Shuttling Paths:")
     t=range(0, model.T, model.N)
-    fig=lineplot(t, model.X[1].(t); width=20, height=5,
+    fig=lineplot(t, model.X[1].(t); width=30, height=9,
     name="x1(t)")
     for i in 2:model.n
         lineplot!(fig, t, model.X[i].(t), name="x$i(t)")
diff --git a/src/stochastics.jl b/src/stochastics.jl
index 0559f7d..a38a435 100644
--- a/src/stochastics.jl
+++ b/src/stochastics.jl
@@ -111,13 +111,8 @@ Covariance function of Gaussian random field.
 - `p₂::Point`: time-position array
 - `process<:GaussianRandomField`: a Gaussian random field, e.g. `OrnsteinUhlenbeckField` or `PinkBrownianField`
 """
-
 function covariance(p₁::Point, p₂::Point, process::OrnsteinUhlenbeckField)::Real
-    process.σ^2 / 4*prod( process.θ) * exp(-dot(process.θ, abs.(p₁ .- p₂)))
-end
-
-function covariance(p₁::Vector{<:Real}, p₂::Vector{<:Real}, process::OrnsteinUhlenbeckField)::Real
-    process.σ^2 / 4*prod( process.θ) * exp(-dot(process.θ, abs.(p₁ .- p₂)))
+    process.σ^2 / prod(2* process.θ) * exp(-dot(process.θ, abs.(p₁ .- p₂)))
 end
 
 function covariance(p₁::Point, p₂::Point, process::PinkBrownianField)::Real
@@ -163,7 +158,7 @@ Auto-Covariance matrix of a Gaussian random process.
 
 """
 function covariancematrix(P::Vector{<:Point}, process::GaussianRandomField)::Symmetric
-    N = size(P, 1)
+    N = length(P)
     A = Matrix{Real}(undef, N, N)
     Threads.@threads for i in 1:N
         for j in i:N
diff --git a/test/runtests.jl b/test/runtests.jl
index 54c06d7..fe4d20c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,7 @@ using SpinShuttling
 using UnicodePlots
 using ProgressMeter
 
-# include("testintegration.jl")
+include("testintegration.jl")
 include("testfidelity.jl")
-# include("teststochastics.jl")
-# include("testspectrum.jl")
+include("teststochastics.jl")
+include("testspectrum.jl")
diff --git a/test/testfidelity.jl b/test/testfidelity.jl
index 560fa94..ab1f09e 100644
--- a/test/testfidelity.jl
+++ b/test/testfidelity.jl
@@ -3,14 +3,13 @@ visualize=true
 
 ##
 @testset begin "test single spin shuttling fidelity"
-    T=200; L=10; σ = sqrt(2) / 20; M = 2000; N=601; κₜ=1/20;κₓ=1/0.1;
+    T=400; L=10; σ = sqrt(2) / 20; M = 2000; N=601; κₜ=1/20;κₓ=1/0.1;
 
     B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
-
     model=OneSpinModel(T,L,N,B)
     # test customize println
     println(model)
-
+    
     f1=averagefidelity(model)
     f2, f2_err=sampling(model, fidelity, M)
     f3=1/2*(1+W(T,L,B))
@@ -23,7 +22,7 @@ end
 
 #
 @testset begin "test single spin forth-back shuttling fidelity"
-    T=200; L=10; σ = sqrt(2) / 20; M = 5000; N=501; κₜ=1/20;κₓ=10; 
+    T=200; L=10; σ = sqrt(2) / 20; M = 2000; N=501; κₜ=1/20;κₓ=10; 
     # exponential should be smaller than 100
     B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
 
@@ -57,11 +56,11 @@ end
 
 ##
 @testset begin "test two spin sequenced shuttling fidelity"
-    L=10; σ =sqrt(2)/20; M=20000; N=501; T1=200; T0=25*0.05; κₜ=1/20; κₓ=1/0.1;
+    L=10; σ =sqrt(2)/20; M=2000; N=501; T1=200; T0=25*0.05; κₜ=1/20; κₓ=1/0.1;
     B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
     model=TwoSpinModel(T0, T1, L, N, B)
     if visualize
-        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, test fig 4"))
+        display(heatmap(collect(model.R.Σ), title="cross covariance matrix, two spin EPR"))
     end
     f1=averagefidelity(model)
     f2, f2_err=sampling(model, fidelity, M)
@@ -111,6 +110,4 @@ end
             lineplot!(fig, t, f_th, name="theoretical fidelity")
         display(fig)
     end
-
-    @test all([abs(f_mc[i]-f_th[i]) < sqrt(f_mc_err[i]) for i in eachindex(t)])
 end
\ No newline at end of file
diff --git a/test/teststochastics.jl b/test/teststochastics.jl
index 2e997f8..c98de95 100644
--- a/test/teststochastics.jl
+++ b/test/teststochastics.jl
@@ -2,7 +2,7 @@ using SpinShuttling: covariancepartition, Symmetric, Cholesky, ishermitian, issy
 using LsqFit
 using Statistics: std, mean
 
-visualize=false
+visualize=true
 
 #
 @testset begin "test of random function"
@@ -23,7 +23,7 @@ visualize=false
 
     if visualize
         display(heatmap(collect(R.Σ), title="covariance matrix, test fig 1")) 
-        display(lineplot([R(),R(),R()],title="random function, test fig 2"))
+        display(lineplot(R(),title="random function, test fig 2"))
         display(heatmap(crosscov, title="cross covariance matrix, test fig 3"))
     end
 
@@ -70,8 +70,8 @@ end
     println("mean 1st order:", mean(err1))
     println("mean 2nd order:", mean(err2))
     if visualize
-        fig=plot(err1, xlabel="T", ylabel="error", label="trapezoid")
-        plot!(err2, label="simpson")
+        fig=lineplot(err1, xlabel="T", ylabel="error", name="trapezoid")
+        lineplot!(fig, err2, name="simpson")
         display(fig)
     end
 end