diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 3fb870af3..8f699f6f5 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -7,6 +7,7 @@ steps:
         - "WeakConvergence2"
         - "WeakConvergence3"
         - "WeakConvergence6"
+        - "WeakConvergence7"
         - "WeakAdaptiveCPU"
     env:
       GROUP: "{{matrix}}"
diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl
index 61cff9659..b58481cbf 100644
--- a/src/StochasticDiffEq.jl
+++ b/src/StochasticDiffEq.jl
@@ -169,7 +169,7 @@ end
 
   export TauLeaping, CaoTauLeaping
 
-  export BAOAB
+  export BAOAB, ABOBA, OBABO
 
   export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm,
          StochasticDiffEqRODECompositeAlgorithm
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index 0d49bd483..1bbb0d036 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -122,6 +122,8 @@ alg_order(alg::TauLeaping) = 1 // 1
 alg_order(alg::CaoTauLeaping) = 1 // 1
 
 alg_order(alg::BAOAB) = 1 // 1
+alg_order(alg::ABOBA) = 2 // 1
+alg_order(alg::OBABO) = 2 // 1
 
 alg_order(alg::SKenCarp) = 2 // 1
 alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs))
@@ -225,6 +227,8 @@ alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::IIF1M) = true
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::IIF2M) = true
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_compatible(prob, a) for a in alg.algs)...)
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::BAOAB) = is_diagonal_noise(prob)
+alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::ABOBA) = is_diagonal_noise(prob)
+alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::OBABO) = is_diagonal_noise(prob)
 
 function alg_compatible(prob::JumpProblem, alg::Union{StochasticDiffEqJumpAdaptiveAlgorithm,StochasticDiffEqJumpAlgorithm})
     prob.prob isa DiscreteProblem
@@ -259,6 +263,7 @@ alg_needs_extra_process(alg::RS2) = true
 alg_needs_extra_process(alg::PL1WM) = true
 alg_needs_extra_process(alg::NON) = true
 alg_needs_extra_process(alg::NON2) = true
+alg_needs_extra_process(alg::OBABO) = true
 
 OrdinaryDiffEq._alg_autodiff(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,FDT,ST,CJ,Controller}) where {CS,AD,FDT,ST,CJ,Controller} = Val{AD}()
 OrdinaryDiffEq._alg_autodiff(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,CJ,Controller}) where {CS,AD,FDT,ST,CJ,Controller} = Val{AD}()
diff --git a/src/algorithms.jl b/src/algorithms.jl
index b38800595..ad5927e29 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -93,7 +93,7 @@ Defaults to solving the Ito problem, but RKMilCommute(interpretation=:Stratonovi
 Uses a 1.5/2.0 error estimate for adaptive time stepping.
 Default: ii_approx=IICommutative() does not approximate the Levy area.
 """
-struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm 
+struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm
   interpretation::Symbol
   ii_approx::T
 end
@@ -869,3 +869,35 @@ struct BAOAB{T} <: StochasticDiffEqAlgorithm
   scale_noise::Bool
 end
 BAOAB(;gamma=1.0, scale_noise=true) = BAOAB(gamma, scale_noise)
+
+
+@doc raw"""
+Leimkuhler B., Matthews C., Robust and efficient configurational molecular sampling via
+Langevin dynamics, J. Chem. Phys. 138, 174102 (2013)
+DOI:10.1063/1.4802990
+
+```math
+du = vdt \\
+dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW
+```
+"""
+struct ABOBA{T} <: StochasticDiffEqAlgorithm
+  gamma::T
+  scale_noise::Bool
+end
+ABOBA(;gamma=1.0, scale_noise=true) = ABOBA(gamma, scale_noise)
+
+@doc raw"""
+Leimkuhler B., Matthews C., Molecular Dynamics With Deterministic and Stochastic Numerical Methods,Interdisciplinary Applied Mathematics; Springer International Publishing: Cham, 2015; Vol. 39
+DOI:10.1007/978-3-319-16375-8
+
+```math
+du = vdt \\
+dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW
+```
+"""
+struct OBABO{T} <: StochasticDiffEqAlgorithm
+  gamma::T
+  scale_noise::Bool
+end
+OBABO(;gamma=1.0, scale_noise=true) = OBABO(gamma, scale_noise)
diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl
index 22dae0cec..895408515 100644
--- a/src/caches/dynamical_caches.jl
+++ b/src/caches/dynamical_caches.jl
@@ -1,42 +1,188 @@
+abstract type StochasticDynamicalEqConstantCache <: StochasticDiffEqConstantCache end # Pourquoi faire ça, Si c'est pour avoir une seul function de check dans initialize!
+abstract type StochasticDynamicalEqMutableCache <: StochasticDiffEqMutableCache end
 
-mutable struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache
+
+mutable struct BAOABConstantCache{uType,uEltypeNoUnits,uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
   k::uType
   half::uEltypeNoUnits
-  c1::uEltypeNoUnits
-  c2::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
 end
-@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType,uTypeCombined} <: StochasticDiffEqMutableCache
+@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
   utmp::uType
+  dumid::uType
   dutmp::uType
+  dunoise::uType
   k::uType
-  gtmp::uType
-  noise::rateNoiseType
+  gtmp::rateNoiseType
+  noise::uType
   half::uEltypeNoUnits
-  c1::uEltypeNoUnits
-  c2::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
   tmp::uTypeCombined
 end
 
 function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
   k = zero(rate_prototype.x[1])
-  c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
-  BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U# if scale_noise == false, c2 = 1
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)# if scale_noise == false, c2 = 1
+  end
+  BAOABConstantCache(k, uEltypeNoUnits(1//2),c2, c2)
 end
 
 function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  dumid = zero(u.x[1])
+  dutmp = zero(u.x[1])
+  dunoise = zero(u.x[1])
+  utmp = zero(u.x[2])
+  k = zero(rate_prototype.x[1])
+
+  gtmp = zero(noise_rate_prototype)
+  noise = zero(rate_prototype.x[1])
+
+  half = uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U# if scale_noise == false, c2 = 1
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)# if scale_noise == false, c2 = 1
+  end
+
+  tmp = zero(u)
+
+  BAOABCache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
+end
+
+
+
+mutable struct ABOBAConstantCache{uType,uEltypeNoUnits, uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
+  k::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+end
+@cache struct ABOBACache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
+  utmp::uType
+  dumid::uType
+  dutmp::uType
+  dunoise::uType
+  k::uType
+  gtmp::rateNoiseType
+  noise::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+  tmp::uTypeCombined
+end
+
+function alg_cache(alg::ABOBA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  k = zero(rate_prototype.x[1])
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+   # if scale_noise == false, c2 = 1
+  ABOBAConstantCache(k, uEltypeNoUnits(1//2), c2, σ)
+end
+
+function alg_cache(alg::ABOBA,prob,u,ΔW,ΔZ,p,rate_prototype, noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
   dutmp = zero(u.x[1])
+  dumid = zero(u.x[1])
+  dunoise = zero(u.x[1])
   utmp = zero(u.x[2])
   k = zero(rate_prototype.x[1])
 
-  gtmp = zero(rate_prototype.x[1])
+  gtmp = zero(noise_rate_prototype)
   noise = zero(rate_prototype.x[1])
 
   half = uEltypeNoUnits(1//2)
-  c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+
+  tmp = zero(u)
+
+  ABOBACache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
+end
+
+
+
+
+mutable struct OBABOConstantCache{uType,rateNoiseType, uEltypeNoUnits, uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
+  k::uType
+  gt::rateNoiseType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+end
+
+@cache struct OBABOCache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
+  utmp::uType
+  dumid::uType
+  dutmp::uType
+  dunoise::uType
+  k::uType
+  gtmp::rateNoiseType
+  noise::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+  tmp::uTypeCombined
+end
+
+function alg_cache(alg::OBABO,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  k = zero(rate_prototype.x[1])
+  gt = zero(noise_rate_prototype)
+  half=uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*half*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*half*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+   # if scale_noise == false, c2 = 1
+  OBABOConstantCache(k, gt, half, c2, σ)
+end
+
+function alg_cache(alg::OBABO,prob,u,ΔW,ΔZ,p,rate_prototype, noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  dutmp = zero(u.x[1])
+  dumid = zero(u.x[1])
+  dunoise = zero(u.x[1])
+  utmp = zero(u.x[2])
+  k = zero(rate_prototype.x[1])
+
+  gtmp = zero(noise_rate_prototype)
+  noise = zero(rate_prototype.x[1])
+
+
+  half = uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*half*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*half*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
 
   tmp = zero(u)
 
-  BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2), tmp)
+  OBABOCache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
 end
diff --git a/src/perform_step/dynamical.jl b/src/perform_step/dynamical.jl
index 074cbc7dd..307bd8ff3 100644
--- a/src/perform_step/dynamical.jl
+++ b/src/perform_step/dynamical.jl
@@ -1,8 +1,8 @@
-function verify_f2(f, p, q, pa, t, integrator, ::BAOABConstantCache)
+function verify_f2(f, p, q, pa, t, integrator, ::StochasticDynamicalEqConstantCache)
     res = f(p, q, pa, t)
     res != p && throwex(integrator)
 end
-function verify_f2(f, res, p, q, pa, t, integrator, ::BAOABCache)
+function verify_f2(f, res, p, q, pa, t, integrator, ::StochasticDiffEqMutableCache)
     f(res, p, q, pa, t)
     res != p && throwex(integrator)
 end
@@ -11,7 +11,7 @@ function throwex(integrator)
   throw(ArgumentError("Algorithm $algn is not applicable if f2(p, q, t) != p"))
 end
 
-function initialize!(integrator, cache::BAOABConstantCache)
+function initialize!(integrator, cache::Union{BAOABConstantCache,ABOBAConstantCache})
   @unpack t,dt,uprev,u,p,W = integrator
   du1 = integrator.uprev.x[1]
   u1 = integrator.uprev.x[2]
@@ -20,7 +20,7 @@ function initialize!(integrator, cache::BAOABConstantCache)
   cache.k = integrator.f.f1(du1,u1,p,t)
 end
 
-function initialize!(integrator, cache::BAOABCache)
+function initialize!(integrator, cache::Union{BAOABCache,ABOBACache})
   @unpack t,dt,uprev,u,p,W = integrator
   du1 = integrator.uprev.x[1]
   u1 = integrator.uprev.x[2]
@@ -31,7 +31,7 @@ end
 
 @muladd function perform_step!(integrator,cache::BAOABConstantCache)
   @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
-  @unpack half, c1, c2 = cache
+  @unpack half, c2, σ = cache
   du1 = uprev.x[1]
   u1 = uprev.x[2]
 
@@ -42,8 +42,12 @@ end
   u2 = u1 + half*dt*du2
 
   # O
-  noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt
-  du3 = c1*du2 + c2*noise
+  noise = integrator.g(u2,p,t+dt*half).*W.dW ./ sqdt
+  if typeof(σ) <: AbstractMatrix || typeof(noise) <: Number
+    du3 = c2*du2 + σ*noise
+  else
+    du3 = c2.*du2 + σ.*noise
+  end
 
   # A
   u = u2 + half*dt*du3
@@ -57,20 +61,26 @@ end
 
 @muladd function perform_step!(integrator,cache::BAOABCache)
   @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
-  @unpack utmp, dutmp, k, gtmp, noise, half, c1, c2 = cache
+  @unpack utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ = cache
   du1 = uprev.x[1]
   u1 = uprev.x[2]
 
   # B
-  @.. dutmp = du1 + half*dt*k
+  @.. dumid = du1 + half*dt*k
 
   # A
-  @.. utmp = u1 + half*dt*dutmp
+  @.. utmp = u1 + half*dt*dumid
 
   # O
   integrator.g(gtmp,utmp,p,t+dt*half)
   @.. noise = gtmp*W.dW / sqdt
-  @.. dutmp = c1*dutmp + c2*noise
+  if typeof(σ) <: AbstractMatrix
+      mul!(dutmp,c2,dumid)
+      mul!(dunoise,σ,noise)
+      @.. dutmp+= dunoise
+  else
+      @.. dutmp = c2*dumid + σ*noise
+  end
 
   # A
   @.. u.x[2] = utmp + half*dt*dutmp
@@ -79,3 +89,168 @@ end
   f.f1(k,dutmp,u.x[2],p,t+dt)
   @.. u.x[1] = dutmp + half*dt*k
 end
+
+
+@muladd function perform_step!(integrator,cache::ABOBAConstantCache)
+  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
+  @unpack half, c2, σ = cache
+  du1 = uprev.x[1]
+  u1 = uprev.x[2]
+
+  # A
+  u_mid = u1 + half*dt*du1
+
+  # BOB: du_t+1
+  cache.k = f.f1(du1,u_mid,p,t+half*dt)
+  noise = integrator.g(u_mid,p,t+dt*half).*W.dW / sqdt
+
+  if typeof(σ) <: AbstractMatrix || typeof(noise) <: Number
+    du = c2 * (du1 + half*dt .* cache.k) .+ σ*noise .+ half * dt .*cache.k
+  else
+    du = c2 .* (du1 + half*dt .* cache.k) .+ σ.*noise .+ half * dt .*cache.k
+  end
+  # A: xt+1
+  u = u_mid .+ half * dt .*du
+
+  integrator.u = ArrayPartition((du, u))
+end
+
+
+@muladd function perform_step!(integrator,cache::ABOBACache)
+  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
+  @unpack utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ = cache
+  du1 = uprev.x[1]
+  u1 = uprev.x[2]
+
+  # A: xt+1/2
+  @.. utmp = u1 + half*dt*du1
+
+
+  # B
+  f.f1(k,du1,utmp,p,t+dt)
+  @.. dumid = du1 + half*dt*k
+
+  # O
+  integrator.g(gtmp,utmp,p,t+dt*half)
+  @.. noise = gtmp*W.dW / sqdt
+
+  if typeof(σ) <: AbstractMatrix
+      mul!(dutmp,c2,dumid)
+      mul!(dunoise,σ,noise)
+      @.. dutmp+=dunoise
+  else
+      @.. dutmp = c2*dumid + σ*noise
+  end
+
+
+  # B
+  @.. u.x[1] = dutmp + half*dt*k
+
+  # A: xt+1
+  @.. u.x[2] = utmp + half*dt*u.x[1]
+end
+
+
+
+
+function initialize!(integrator, cache::OBABOConstantCache)
+  @unpack t,dt,uprev,u,p,W = integrator
+  du1 = integrator.uprev.x[1]
+  u1 = integrator.uprev.x[2]
+
+  verify_f2(integrator.f.f2, du1, u1, p, t, integrator, cache)
+  cache.k = integrator.f.f1(du1,u1,p,t)
+  cache.gt = integrator.g(u1,p,t)
+end
+
+function initialize!(integrator, cache::OBABOCache)
+  @unpack t,dt,uprev,u,p,W = integrator
+  du1 = integrator.uprev.x[1]
+  u1 = integrator.uprev.x[2]
+
+  verify_f2(integrator.f.f2, cache.k, du1, u1, p, t, integrator, cache)
+  integrator.f.f1(cache.k,du1,u1,p,t)
+  integrator.g(cache.gtmp,u1,p,t)
+end
+
+
+@muladd function perform_step!(integrator,cache::OBABOConstantCache)
+  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
+  @unpack half, c2, σ = cache
+  du1 = uprev.x[1]
+  u1 = uprev.x[2]
+
+  # O
+  noise = cache.gt.*W.dW ./ sqdt
+  if typeof(σ) <: AbstractMatrix || typeof(noise) <: Number
+    du2 = c2*du1 + σ*noise
+  else
+    du2 = c2.*du1 + σ.*noise
+  end
+
+  # B
+  dumid = du2 + half*dt*cache.k
+
+  # A
+  u = u1 + dt*dumid
+
+  cache.k = f.f1(dumid,u,p,t+dt)
+  # B
+  du3 = dumid + half*dt*cache.k
+
+  # O
+  cache.gt = integrator.g(u,p,t+dt)
+  noise = cache.gt.*W.dZ ./ sqdt # That should be a second noise
+  if typeof(σ) <: AbstractMatrix || typeof(noise) <: Number
+    du = c2*du3 + σ*noise
+  else
+    du = c2.*du3 + σ.*noise
+  end
+
+  integrator.u = ArrayPartition((du, u))
+end
+
+
+@muladd function perform_step!(integrator,cache::OBABOCache)
+  @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
+  @unpack utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ = cache
+  du1 = uprev.x[1]
+  u1 = uprev.x[2]
+
+  # O
+  @.. noise = gtmp*W.dW / sqdt
+
+  if typeof(σ) <: AbstractMatrix
+      mul!(dutmp,c2,du1)
+      mul!(dunoise,σ,noise)
+      @.. dutmp+=dunoise
+  else
+      @.. dutmp = c2*du1 + σ*noise
+  end
+
+  # B
+
+  @.. dumid = dutmp + half*dt*k
+
+  # A: xt+1
+  @.. u.x[2] = u1 + dt*dumid
+
+
+  # B
+  f.f1(k,dumid,u.x[2],p,t+dt)
+  @..  dutmp = dumid + half*dt*k
+
+  # O
+  integrator.g(gtmp,u.x[2],p,t+dt)
+  @.. noise = gtmp*W.dZ / sqdt  # That should be a second noise
+
+  if typeof(σ) <: AbstractMatrix
+      mul!(u.x[1],c2,dutmp)
+      mul!(dunoise,σ,noise)
+      @.. u.x[1]+=dunoise
+  else
+      @.. u.x[1] = c2*dutmp + σ*noise
+  end
+
+
+end
diff --git a/src/solve.jl b/src/solve.jl
index c387fb16a..daaa4b262 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -211,10 +211,13 @@ function DiffEqBase.__init(
   end
   rateType = typeof(rate_prototype) ## Can be different if united
 
-  if prob.f isa DynamicalSDEFunction
-    noise_rate_prototype = rate_prototype.x[1]
-  elseif is_diagonal_noise(prob)
-    noise_rate_prototype = rate_prototype
+
+  if is_diagonal_noise(prob)
+    if prob.f isa DynamicalSDEFunction
+      noise_rate_prototype = rate_prototype.x[1]
+    else
+      noise_rate_prototype = rate_prototype
+    end
   elseif prob isa DiffEqBase.AbstractRODEProblem
     if prob isa DiffEqBase.AbstractSDEProblem
       noise_rate_prototype = copy(prob.noise_rate_prototype)
@@ -288,7 +291,7 @@ function DiffEqBase.__init(
   else
     randElType = uBottomEltypeNoUnits # Strip units and type info
     if prob.f isa DynamicalSDEFunction
-      rand_prototype = copy(noise_rate_prototype)
+      rand_prototype = copy(rate_prototype.x[1]) # Noise is a vector of the same size than the number of variables
     elseif is_diagonal_noise(prob)
       if u isa SArray
         rand_prototype = zero(u) # TODO: Array{randElType} for units
diff --git a/test/runtests.jl b/test/runtests.jl
index abe68baa4..0dd7aa910 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -41,6 +41,7 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
     @time @safetestset "Stiffness Detection Test" begin include("stiffness_detection_test.jl") end
     @time @safetestset "Adaptive SDE Linear Tests" begin include("adaptive/sde_linearadaptive_tests.jl") end
     @time @safetestset "Inplace RODESolution Interpolation Tests" begin include("inplace_interpolation.jl") end
+    @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end
   end
 
   if GROUP == "All" || GROUP == "Interface3"
@@ -59,7 +60,6 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
 
   if !is_APPVEYOR && (GROUP == "All" || GROUP == "AlgConvergence")
     @time @safetestset "Convergence Tests" begin include("sde/sde_convergence_tests.jl") end
-    @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end
   end
 
   if !is_APPVEYOR && GROUP == "AlgConvergence2"
@@ -102,6 +102,10 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
     @time @safetestset "Roessler weak SRK diagonal Tests" begin include("weak_convergence/srk_weak_diagonal_final.jl") end
   end
 
+  if !is_APPVEYOR && GROUP == "WeakConvergence7"
+     @time @safetestset "Dynamical SDE Tests" begin include("weak_convergence/sde_convergence_dynamical.jl") end
+   end
+
   if !is_APPVEYOR && GROUP == "OOPWeakConvergence"
     @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end
     @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end
diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl
index a0a259810..d94b13baa 100644
--- a/test/sde/sde_dynamical.jl
+++ b/test/sde/sde_dynamical.jl
@@ -3,45 +3,77 @@ Random.seed!(1)
 
 f1_harmonic(v,u,p,t) = -u
 f2_harmonic(v,u,p,t) = v
-g(u,p,t) = 1
+g(u,p,t) = 1 .+zero(u)
 γ = 1
 
-@testset "Scalar u" begin
-    u0 = 0
-    v0 = 1
+
+@testset "Vector u" begin
+
+    u0 = zeros(2)
+    v0 = ones(2)
+
+    f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
+    f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
+    g_iip(du,u,p,t) = du .= g(u,p,t)
 
     ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
     prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
+    sol1 = solve(prob1,BAOAB(gamma=[γ,γ]);dt=1/10,save_noise=true)
+
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,BAOAB(gamma=[γ,γ]);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
 
     dts = (1/2) .^ (8:-1:4)
 
-    # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
-    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(2e2),use_noise_grid=false)
-    display(sim1.𝒪est)
-    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+    sol1 = solve(prob1,ABOBA(gamma=[γ,γ]);dt=1/10,save_noise=true)
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,ABOBA(gamma=[γ,γ]);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
+
+    sol1 = solve(prob1,OBABO(gamma=[γ,γ]);dt=1/10,save_noise=true)
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,OBABO(gamma=[γ,γ]);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
+
 end
 
-@testset "Vector u" begin
+
+
+@testset "Matricial gamma" begin
 
     u0 = zeros(2)
     v0 = ones(2)
 
+    gamma_mat = [γ -0.1*γ;0.5*γ γ]
+
     f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
     f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
     g_iip(du,u,p,t) = du .= g(u,p,t)
 
     ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
     prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
-    sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true)
+    sol1 = solve(prob1,BAOAB(gamma=gamma_mat);dt=1/10,save_noise=true)
 
     prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
-    sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10)
+    sol2 = solve(prob2,BAOAB(gamma=gamma_mat);dt=1/10)
 
     @test sol1[:] ≈ sol2[:]
 
-    dts = (1/2) .^ (8:-1:4)
 
-    # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
-    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
-    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+    sol1 = solve(prob1,ABOBA(gamma=gamma_mat);dt=1/10,save_noise=true)
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,ABOBA(gamma=gamma_mat);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
+
+
+    sol1 = solve(prob1,OBABO(gamma=gamma_mat);dt=1/10,save_noise=true)
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,OBABO(gamma=gamma_mat);dt=1/10)
+
+    @test sol1[:] ≈ sol2[:]
 end
diff --git a/test/weak_convergence/sde_convergence_dynamical.jl b/test/weak_convergence/sde_convergence_dynamical.jl
new file mode 100644
index 000000000..5d604fea7
--- /dev/null
+++ b/test/weak_convergence/sde_convergence_dynamical.jl
@@ -0,0 +1,31 @@
+using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random
+Random.seed!(1)
+
+f1_harmonic(v,u,p,t) = -u
+f2_harmonic(v,u,p,t) = v
+g(u,p,t) = 1 .+zero(u)
+γ = 1
+
+@testset "Dynamical convergence" begin
+    u0 = 0
+    v0 = 1
+
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,1.5))
+
+    dts = (1/2) .^ (6:-1:3)
+
+    # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
+    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e5),use_noise_grid=false,verbose=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-2) < 0.5
+
+    sim1  = analyticless_test_convergence(dts,prob1,ABOBA(gamma=γ),(1/2)^10;trajectories=Int(1e5),use_noise_grid=false,verbose=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-2) < 0.5
+
+
+    sim1  = analyticless_test_convergence(dts,prob1,OBABO(gamma=γ),(1/2)^10;trajectories=Int(1e5),use_noise_grid=false,verbose=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-2) < 0.5
+end