diff --git a/Project.toml b/Project.toml
index c5a003c4b3..d731092df7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.6.9"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
-FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
 Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
@@ -29,7 +28,6 @@ Colors = "0.12"
 Distributions = "0.22.6, 0.23, 0.24, 0.25"
 Einsum = "0.4"
 FiniteDiff = "2"
-FiniteDifferences = "0.12"
 HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 LightGraphs = "1"
@@ -48,6 +46,7 @@ julia = "1.5"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
@@ -65,4 +64,4 @@ VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "Zygote"]
+test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "FiniteDifferences", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "Zygote"]
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 55b5b817d4..f89f3246cc 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -86,7 +86,6 @@ import Base:
 using Base.Iterators: repeated
 using Distributions
 using Einsum: @einsum
-using FiniteDifferences
 using HybridArrays
 using Kronecker
 using LightGraphs
@@ -287,6 +286,11 @@ function __init__()
         include("differentiation/finite_diff.jl")
     end
 
+    @require FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" begin
+        using .FiniteDifferences
+        include("differentiation/finite_differences.jl")
+    end
+
     @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
         using .ForwardDiff
         include("differentiation/forward_diff.jl")
@@ -294,7 +298,7 @@ function __init__()
 
     @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
         using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
-        include("ode.jl")
+        include("differentiation/ode.jl")
     end
 
     @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
@@ -440,6 +444,7 @@ export AbstractRetractionMethod,
     PolarRetraction,
     ProjectionRetraction,
     SoftmaxRetraction,
+    ODEExponentialRetraction,
     PadeRetraction,
     ProductRetraction,
     PowerRetraction
@@ -664,9 +669,9 @@ export get_basis,
 export AbstractDiffBackend,
     AbstractRiemannianDiffBackend,
     FiniteDifferencesBackend,
-    RiemannianONBDiffBackend,
-    RiemannianProjectionGradientBackend
-export diff_backend, diff_backend!, diff_backends
+    TangentDiffBackend,
+    RiemannianProjectionBackend
+export default_differential_backend, set_default_differential_backend!
 # atlases and charts
 export get_point, get_point!, get_parameters, get_parameters!
 
diff --git a/src/differentiation/differentiation.jl b/src/differentiation/differentiation.jl
index 17ed083186..45d12c177a 100644
--- a/src/differentiation/differentiation.jl
+++ b/src/differentiation/differentiation.jl
@@ -14,7 +14,7 @@ struct NoneDiffBackend <: AbstractDiffBackend end
 
 Compute the derivative of a callable `f` at time `t` computed using the given `backend`,
 an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly
-specified, it is obtained using the function [`Manifolds.diff_backend`](@ref).
+specified, it is obtained using the function [`default_differential_backend`](@ref).
 
 This function calculates plain Euclidean derivatives, for Riemannian differentiation see
 for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)).
@@ -26,7 +26,9 @@ for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::A
 """
 function _derivative end
 
-function _derivative!(f, X, t, backend::AbstractDiffBackend)
+_derivative(f, t) = _derivative(f, t, default_differential_backend())
+
+function _derivative!(f, X, t, backend::AbstractDiffBackend=default_differential_backend())
     return copyto!(X, _derivative(f, t, backend))
 end
 
@@ -35,7 +37,7 @@ end
 
 Compute the gradient of a callable `f` at point `p` computed using the given `backend`,
 an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
-specified, it is obtained using the function [`diff_backend`](@ref).
+specified, it is obtained using the function [`default_differential_backend`](@ref).
 
 This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
 for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).
@@ -47,98 +49,132 @@ for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::An
 """
 function _gradient end
 
-function _gradient!(f, X, p, backend::AbstractDiffBackend)
+_gradient(f, p) = _gradient(f, p, default_differential_backend())
+
+function _gradient!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
     return copyto!(X, _gradient(f, p, backend))
 end
 
+"""
+    _jacobian(f, p[, backend::AbstractDiffBackend])
+
+Compute the jacobian of a callable `f` at point `p` computed using the given `backend`,
+an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
+specified, it is obtained using the function [`default_differential_backend`](@ref).
+
+This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
+for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).
+
+!!! note
+
+    Not specifying the backend explicitly will usually result in a type instability
+    and decreased performance.
+"""
+function _jacobian end
+
+_jacobian(f, p) = _jacobian(f, p, default_differential_backend())
+
+function _jacobian!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
+    return copyto!(X, _jacobian(f, p, backend))
+end
+
 """
     CurrentDiffBackend(backend::AbstractDiffBackend)
 
 A mutable struct for storing the current differentiation backend in a global
-constant [`_current_diff_backend`](@ref).
+constant [`_current_default_differential_backend`](@ref).
 
 # See also
 
-[`AbstractDiffBackend`](@ref), [`diff_backend`](@ref), [`diff_backend!`](@ref)
+[`AbstractDiffBackend`](@ref), [`default_differential_backend`](@ref), [`set_default_differential_backend`](@ref)
 """
 mutable struct CurrentDiffBackend
     backend::AbstractDiffBackend
 end
 
 """
-    _current_diff_backend
+    _current_default_differential_backend
 
 The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default
 differentiation backend.
 """
-const _current_diff_backend = CurrentDiffBackend(NoneDiffBackend())
-
+const _current_default_differential_backend = CurrentDiffBackend(NoneDiffBackend())
 """
-    _diff_backends
+    default_differential_backend() -> AbstractDiffBackend
 
-A vector of valid [`Manifolds.AbstractDiffBackend`](@ref).
+Get the default differentiation backend.
 """
-const _diff_backends = AbstractDiffBackend[]
+default_differential_backend() = _current_default_differential_backend.backend
 
 """
-    diff_backend() -> AbstractDiffBackend
-
-Get the current differentiation backend.
-"""
-diff_backend() = _current_diff_backend.backend
-
-"""
-    diff_backend!(backend::AbstractDiffBackend)
+    set_default_differential_backend!(backend::AbstractDiffBackend)
 
 Set current backend for differentiation to `backend`.
 """
-function diff_backend!(backend::AbstractDiffBackend)
-    _current_diff_backend.backend = backend
+function set_default_differential_backend!(backend::AbstractDiffBackend)
+    _current_default_differential_backend.backend = backend
     return backend
 end
 
-"""
-    diff_backends() -> Vector{AbstractDiffBackend}
-
-Get vector of currently valid differentiation backends.
-"""
-diff_backends() = _diff_backends
-
-_derivative(f, t) = _derivative(f, t, diff_backend())
-
-_derivative!(f, X, t) = _derivative!(f, X, t, diff_backend())
+@doc raw"""
+    ODEExponentialRetraction{T<:AbstractRetractionMethod, B<: AbstractBasis} <: AbstractRetractionMethod
 
-_gradient(f, p) = _gradient(f, p, diff_backend())
+This retraction approximates the exponential map by solving the correspondig ODE.
+Let ``p\in\mathal M`` and ``X\in T_p\mathcal M`` denote the input for the exponential map
+and ``d`` denote the [`manifold_dimension`](@ref) of `M`.
 
-_gradient!(f, X, p) = _gradient!(f, X, p, diff_backend())
+This the ODE is formulated in a chart constructed using an [`AbstractBasis`](@ref) `B` and an
+[`AbstractRetractionMethod`](@ref) `R` as follows.
+Given some coordinates ``c\in ℝ^d`` - these can be used to form a tangent vector
+with restect to th basis `B, i.e. ``c \mapsto Y=``[`get_vector`](@ref)`(M, p, c, B)`.
+Further, using the retraction we can map ``Y`` to a point on the manifold
+``Y \mapsto q =``[`retract`](@ref)`(M, p, X, R)`.
 
-# Finite differences
-
-"""
-    FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1))
-
-Differentiation backend based on the FiniteDifferences package.
+Hence the ODE can be formulated in a curve ``c(t)`` in parameter space.
+This is – for sure – only possible locally as fas as the retraction is well-defined.
 """
-struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend
-    method::TM
+struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <:
+       AbstractRetractionMethod
+    retraction::T
+    basis::B
 end
-
-function FiniteDifferencesBackend()
-    return FiniteDifferencesBackend(central_fdm(5, 1))
+function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod}
+    return ODEExponentialRetraction(r, DefaultOrthonormalBasis())
 end
-
-push!(_diff_backends, FiniteDifferencesBackend())
-
-diff_backend!(_diff_backends[end])
-
-function _derivative(f, t, backend::FiniteDifferencesBackend)
-    return backend.method(f, t)
+function ODEExponentialRetraction(r::T, ::CachedBasis) where {T<:AbstractRetractionMethod}
+    return throw(
+        DomainError(
+            r,
+            "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.",
+        ),
+    )
 end
-
-function _gradient(f, p, backend::FiniteDifferencesBackend)
-    return FiniteDifferences.grad(backend.method, f, p)[1]
+function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis)
+    return throw(
+        DomainError(
+            r,
+            "You can not use the exponential map as an inner method to solve the ode for the exponential map.",
+        ),
+    )
+end
+function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis)
+    return throw(
+        DomainError(
+            r,
+            "Neither the exponential map nor a Cached Basis can be used with this retraction type.",
+        ),
+    )
 end
 
-function _jacobian(f, p, backend::FiniteDifferencesBackend)
-    return FiniteDifferences.jacobian(backend.method, f, p)[1]
+function retract!(M::AbstractManifold, q, p, X, r::ODEExponentialRetraction)
+    sol = solve_exp_ode(
+        M,
+        p,
+        X;
+        basis=r.basis,
+        retraction=r.retraction,
+        dense=false,
+        saveat=[1.0],
+    )
+    return copyto!(q, sol)
 end
diff --git a/src/differentiation/finite_diff.jl b/src/differentiation/finite_diff.jl
index 3599522a4d..bd5403d60a 100644
--- a/src/differentiation/finite_diff.jl
+++ b/src/differentiation/finite_diff.jl
@@ -1,8 +1,11 @@
 
 """
-    FiniteDiffBackend(method::Val{Symbol} = Val{:central})
+    FiniteDiffBackend <: AbstractDiffBackend
+
+A type to specify / use differentiation backend based on FiniteDiff package.
 
-Differentiation backend based on FiniteDiff package.
+# Constructor
+    FiniteDiffBackend(method::Val{Symbol} = Val{:central})
 """
 struct FiniteDiffBackend{TM<:Val} <: AbstractDiffBackend
     method::TM
@@ -10,16 +13,26 @@ end
 
 FiniteDiffBackend() = FiniteDiffBackend(Val(:central))
 
-push!(_diff_backends, FiniteDiffBackend())
-
-function _derivative(f, p, backend::FiniteDiffBackend{Method}) where {Method}
+function _derivative(f, p, ::FiniteDiffBackend{Method}) where {Method}
     return FiniteDiff.finite_difference_derivative(f, p, Method)
 end
 
-function _gradient(f, p, backend::FiniteDiffBackend{Method}) where {Method}
+function _gradient(f, p, ::FiniteDiffBackend{Method}) where {Method}
     return FiniteDiff.finite_difference_gradient(f, p, Method)
 end
 
-function _gradient!(f, X, p, backend::FiniteDiffBackend{Method}) where {Method}
+function _gradient!(f, X, p, ::FiniteDiffBackend{Method}) where {Method}
     return FiniteDiff.finite_difference_gradient!(X, f, p, Method)
 end
+
+function _jacobian(f, p, ::FiniteDiffBackend{Method}) where {Method}
+    return FiniteDiff.finite_difference_jacobian(f, p, Method)
+end
+
+function _jacobian!(f, X, p, ::FiniteDiffBackend{Method}) where {Method}
+    return FiniteDiff.finite_difference_jacobian!(X, f, p, Method)
+end
+
+if default_differential_backend() === NoneDiffBackend()
+    set_default_differential_backend!(FiniteDiffBackend())
+end
diff --git a/src/differentiation/finite_differences.jl b/src/differentiation/finite_differences.jl
new file mode 100644
index 0000000000..6129218c28
--- /dev/null
+++ b/src/differentiation/finite_differences.jl
@@ -0,0 +1,28 @@
+"""
+    FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1))
+
+Differentiation backend based on the FiniteDifferences package.
+"""
+struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend
+    method::TM
+end
+
+function FiniteDifferencesBackend()
+    return FiniteDifferencesBackend(central_fdm(5, 1))
+end
+
+function _derivative(f, t, backend::FiniteDifferencesBackend)
+    return backend.method(f, t)
+end
+
+function _gradient(f, p, backend::FiniteDifferencesBackend)
+    return FiniteDifferences.grad(backend.method, f, p)[1]
+end
+
+function _jacobian(f, p, backend::FiniteDifferencesBackend)
+    return FiniteDifferences.jacobian(backend.method, f, p)[1]
+end
+
+if default_differential_backend() === NoneDiffBackend()
+    set_default_differential_backend!(FiniteDifferencesBackend())
+end
diff --git a/src/differentiation/forward_diff.jl b/src/differentiation/forward_diff.jl
index d6aabe3583..4a06b927ef 100644
--- a/src/differentiation/forward_diff.jl
+++ b/src/differentiation/forward_diff.jl
@@ -21,4 +21,6 @@ function _jacobian(f, p, ::ForwardDiffBackend)
     return ForwardDiff.jacobian(f, p)
 end
 
-push!(_diff_backends, ForwardDiffBackend())
+if default_differential_backend() === NoneDiffBackend()
+    set_default_differential_backend!(ForwardDiffBackend())
+end
diff --git a/src/differentiation/ode.jl b/src/differentiation/ode.jl
new file mode 100644
index 0000000000..4b6dcdb412
--- /dev/null
+++ b/src/differentiation/ode.jl
@@ -0,0 +1,35 @@
+function solve_exp_ode(
+    M::AbstractManifold,
+    p,
+    X;
+    basis::AbstractBasis=DefaultOrthonormalBasis(),
+    solver=AutoVern9(Rodas5()),
+    backend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
+    kwargs...,
+)
+    n = length(p)
+    iv = SVector{n}(1:n)
+    ix = SVector{n}((n + 1):(2 * n))
+    u0 = allocate(p, 2 * n)
+    u0[iv] .= get_get_coordinates!(M, u0[iv])
+    u0[ix] .= p
+
+    function exp_problem(u, params, t)
+        M = params[1]
+        q = u[ix]
+        dx = u[iv]
+        ddx = allocate(u, Size(n))
+        du = allocate(u)
+        Γ = christoffel_symbols_second(M, q, basis; backend=backend, retraction=retraction)
+        @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j]
+        du[iv] .= ddx
+        du[ix] .= dx
+        return Base.convert(typeof(u), du)
+    end
+    params = (M,)
+    prob = ODEProblem(exp_problem, u0, (0.0, 1.0), params)
+    sol = solve(prob, solver; kwargs...)
+    q = sol.u[1][(n + 1):(2 * n)]
+    return q
+end
diff --git a/src/differentiation/reverse_diff.jl b/src/differentiation/reverse_diff.jl
index 851bf89a2d..4ec0081c01 100644
--- a/src/differentiation/reverse_diff.jl
+++ b/src/differentiation/reverse_diff.jl
@@ -8,4 +8,6 @@ function Manifolds._gradient!(f, X, p, ::ReverseDiffBackend)
     return ReverseDiff.gradient!(X, f, p)
 end
 
-push!(Manifolds._diff_backends, ReverseDiffBackend())
+if default_differential_backend() === NoneDiffBackend()
+    set_default_differential_backend!(ReverseDiffBackend())
+end
diff --git a/src/differentiation/riemannian_diff.jl b/src/differentiation/riemannian_diff.jl
index d3cc0d1bea..63395992ae 100644
--- a/src/differentiation/riemannian_diff.jl
+++ b/src/differentiation/riemannian_diff.jl
@@ -2,25 +2,30 @@
 """
     AbstractRiemannianDiffBackend
 
-An abstract type for diff backends. See [`RiemannianONBDiffBackend`](@ref) for
-an example.
+An abstract type for backends for differentiation.
 """
 abstract type AbstractRiemannianDiffBackend end
 
 @doc raw"""
-    differential(M::AbstractManifold, f, t::Real, backend::AbstractDiffBackend = rdifferential_backend())
+    differential(M::AbstractManifold, f, t::Real, backend::AbstractDiffBackend)
+    differential!(M::AbstractManifold, f, X, t::Real, backend::AbstractDiffBackend)
 
 Compute the Riemannian differential of a curve $f: ℝ\to M$ on a manifold `M`
 represented by function `f` at time `t` using the given backend.
 It is calculated as the tangent vector equal to $\mathrm{d}f_t(t)[1]$.
+
+The mutating variant computes the differential in place of `X`.
 """
 differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)
 
 @doc raw"""
-    gradient(M::AbstractManifold, f, p, backend::AbstractRiemannianDiffBackend = rgradient_backend())
+    gradient(M::AbstractManifold, f, p, backend::AbstractRiemannianDiffBackend)
+    gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiffBackend)
+
+Compute the Riemannian gradient ``∇f(p)`` of a real-valued function ``f:\mathcal M \to ℝ``
+at at point `p` on the manifold `M` using the specified [`AbstractRiemannianDiffBackend`](@ref).
 
-Compute the Riemannian gradient $∇f(p)$ of a real field on manifold `M` represented by
-function `f` at point `p` using the given backend.
+The mutating variant computes the gradient in place of `X`.
 """
 gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)
 
@@ -38,43 +43,61 @@ function gradient!(M::AbstractManifold, f, X, p, backend::AbstractRiemannianDiff
     return copyto!(X, gradient(M, f, p, backend))
 end
 
-differential(M::AbstractManifold, f, p) = differential(M, f, p, rdifferential_backend())
+@doc raw"""
+    TangentDiffBackend <: AbstractRiemannianDiffBackend
 
-function differential!(M::AbstractManifold, f, X, p)
-    return differential!(M, f, X, p, rdifferential_backend())
-end
+A backend that uses a tangent space and a basis therein to derive an
+(intrinsic, approximate) differentiation scheme.
 
-gradient(M::AbstractManifold, f, p) = gradient(M, f, p, rgradient_backend())
+Since it works in a tangent space, methods might require a retraction and an
+inverse retraction as well as a basis.
 
-gradient!(M::AbstractManifold, f, X, p) = gradient!(M, f, X, p, rgradient_backend())
+In the tangent space itself, this backend then employs an (Euclidean)
+[`AbstractDiffBackend`](@ref)
 
+# Constructor
+
+    TangentDiffBackend(diff_backend)
+
+where `diff_backend` is an [`AbstractDiffBackend`](@ref) to be used on the tangent space.
+
+With the keyword arguments
+
+* `retraction` an [`AbstractRetractionMethod`](@ref) ([`ExponentialRetraction`]('ref) by default)
+* `inverse_retraction` an [`AbstractInverseRetractionMethod`](@ref) ([`LogarithmicInverseRetraction`]('ref) by default)
+* `basis` an [`AbstractBasis`](@ref) ([`DefaultOrthogonalBasis`]('ref) by default)
 """
-    RiemannianONBDiffBackend(
-        diff_backend::AbstractDiffBackend
-        retraction::AbstractRetractionMethod
-        inverse_retraction::AbstractInverseRetractionMethod
-        basis::Union{AbstractOrthonormalBasis,CachedBasis{<:AbstractOrthonormalBasis}},
-    ) <: AbstractRiemannianDiffBackend
-
-Riemannian differentiation based on differentiation in an [`AbstractOrthonormalBasis`](@ref)
-`basis` using specified `retraction`, `inverse_retraction` and using backend `diff_backend`.
-"""
-struct RiemannianONBDiffBackend{
-    TADBackend<:AbstractDiffBackend,
-    TRetr<:AbstractRetractionMethod,
-    TInvRetr<:AbstractInverseRetractionMethod,
-    TBasis<:Union{
-        AbstractOrthonormalBasis,
-        CachedBasis{𝔽,<:AbstractOrthonormalBasis{𝔽}} where {𝔽},
-    },
+struct TangentDiffBackend{
+    TAD<:AbstractDiffBackend,
+    TR<:AbstractRetractionMethod,
+    TIR<:AbstractInverseRetractionMethod,
+    TB<:AbstractBasis,
 } <: AbstractRiemannianDiffBackend
-    diff_backend::TADBackend
-    retraction::TRetr
-    inverse_retraction::TInvRetr
-    basis::TBasis
+    diff_backend::TAD
+    retraction::TR
+    inverse_retraction::TIR
+    basis::TB
+end
+function TangentDiffBackend(
+    diff_backend::TAD;
+    retraction::TR=ExponentialRetraction(),
+    inverse_retraction::TIR=LogarithmicInverseRetraction(),
+    basis::TB=DefaultOrthonormalBasis(),
+) where {
+    TAD<:AbstractDiffBackend,
+    TR<:AbstractRetractionMethod,
+    TIR<:AbstractInverseRetractionMethod,
+    TB<:AbstractBasis,
+}
+    return TangentDiffBackend{TAD,TR,TIR,TB}(
+        diff_backend,
+        retraction,
+        inverse_retraction,
+        basis,
+    )
 end
 
-function differential(M::AbstractManifold, f, t::Real, backend::RiemannianONBDiffBackend)
+function differential(M::AbstractManifold, f, t::Real, backend::TangentDiffBackend)
     p = f(t)
     onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
         return get_coordinates(
@@ -87,13 +110,7 @@ function differential(M::AbstractManifold, f, t::Real, backend::RiemannianONBDif
     return get_vector(M, p, onb_coords, backend.basis)
 end
 
-function differential!(
-    M::AbstractManifold,
-    f,
-    X,
-    t::Real,
-    backend::RiemannianONBDiffBackend,
-)
+function differential!(M::AbstractManifold, f, X, t::Real, backend::TangentDiffBackend)
     p = f(t)
     onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
         return get_coordinates(
@@ -106,7 +123,35 @@ function differential!(
     return get_vector!(M, X, p, onb_coords, backend.basis)
 end
 
-function gradient(M::AbstractManifold, f, p, backend::RiemannianONBDiffBackend)
+@doc raw"""
+    gradient(M, f, p, backend::TangentDiffBackend)
+
+This method uses the internal `backend.diff_backend` (Euclidean) on the function
+
+```math
+    f(\retr_p(\cdot))
+```
+
+which is given on the tangent space. In detail, the gradient can be written in
+terms of the `backend.basis`. We illustrate it here for an [`AbstractOrthonormalBasis`](@ref),
+since that simplifies notations:
+
+```math
+\operatorname{grad}f(p) = \operatorname{grad}f(p) = \sum_{i=1}^{d} g_p(\operatorname{grad}f(p),X_i)X_i
+	= \sum_{i=1}^{d} Df(p)[X_i]X_i
+```
+
+where the last equality is due to the definition of the gradient as the Riesz representer of the differential.
+
+If the backend is a forward (or backward) finite difference, these coefficients in the sum can be approximates as
+
+```math
+DF(p)[Y] ≈ \frac{1}{h}\bigl( f(\exp_p(hY)) - f(p) \bigr)
+```
+writing ``p=\exp_p(0)`` we see that this is a finite difference of ``f\circ\exp_p``, i.e. for
+a function on the tangent space, so we can also use other (Euclidean) backends
+"""
+function gradient(M::AbstractManifold, f, p, backend::TangentDiffBackend)
     X = get_coordinates(M, p, zero_vector(M, p), backend.basis)
     onb_coords = _gradient(X, backend.diff_backend) do Y
         return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
@@ -114,7 +159,7 @@ function gradient(M::AbstractManifold, f, p, backend::RiemannianONBDiffBackend)
     return get_vector(M, p, onb_coords, backend.basis)
 end
 
-function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianONBDiffBackend)
+function gradient!(M::AbstractManifold, f, X, p, backend::TangentDiffBackend)
     X2 = get_coordinates(M, p, zero_vector(M, p), backend.basis)
     onb_coords = _gradient(X2, backend.diff_backend) do Y
         return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
@@ -122,126 +167,44 @@ function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianONBDiffBacke
     return get_vector!(M, X, p, onb_coords, backend.basis)
 end
 
-"""
-    CurrentRiemannianDiffBackend(backend::AbstractRiemannianDiffBackend)
-
-A mutable struct for storing the current Riemannian differentiation backend in global
-constants [`Manifolds._current_rgradient_backend`](@ref) and
-[`Manifolds._current_rdifferential_backend`](@ref).
-
-# See also
-
-[`AbstractRiemannianDiffBackend`](@ref), [`rdifferential_backend`](@ref),
-[`rdifferential_backend!`](@ref)
-"""
-mutable struct CurrentRiemannianDiffBackend
-    backend::AbstractRiemannianDiffBackend
-end
-
-"""
-    _current_rgradient_backend
-
-The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the
-globally default differentiation backend for calculating gradients.
-
-# See also
-
-[`Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)`](@ref)
-"""
-const _current_rgradient_backend = CurrentRiemannianDiffBackend(
-    RiemannianONBDiffBackend(
-        diff_backend(),
-        ExponentialRetraction(),
-        LogarithmicInverseRetraction(),
-        DefaultOrthonormalBasis(),
-    ),
-)
-
-"""
-    _current_rdifferential_backend
-
-The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the
-globally default differentiation backend for calculating differentials.
-
-# See also
-
-[`Manifolds.differential`](@ref)
-"""
-const _current_rdifferential_backend = CurrentRiemannianDiffBackend(
-    RiemannianONBDiffBackend(
-        diff_backend(),
-        ExponentialRetraction(),
-        LogarithmicInverseRetraction(),
-        DefaultOrthonormalBasis(),
-    ),
-)
-
-"""
-    rgradient_backend() -> AbstractRiemannianDiffBackend
-
-Get the current differentiation backend for Riemannian gradients.
-"""
-rgradient_backend() = _current_rgradient_backend.backend
-
-"""
-    rgradient_backend!(backend::AbstractRiemannianDiffBackend)
-
-Set current Riemannian gradient backend for differentiation to `backend`.
-"""
-function rgradient_backend!(backend::AbstractRiemannianDiffBackend)
-    _current_rgradient_backend.backend = backend
-    return backend
-end
-
-"""
-    rdifferential_backend() -> AbstractRiemannianDiffBackend
-
-Get the current differentiation backend for Riemannian differentials.
-"""
-rdifferential_backend() = _current_rdifferential_backend.backend
-
-"""
-    rdifferential_backend!(backend::AbstractRiemannianDiffBackend)
-
-Set current Riemannian differential backend for differentiation to `backend`.
-"""
-function rdifferential_backend!(backend::AbstractRiemannianDiffBackend)
-    _current_rdifferential_backend.backend = backend
-    return backend
-end
-
-"""
-    RiemannianProjectionGradientBackend(
-        diff_backend::AbstractDiffBackend
-    ) <: AbstractRiemannianDiffBackend
-
-Riemannian differentiation based on differentiation in the ambient space and projection to
-the given manifold. Differentiation in the ambient space is performed using
-the backend `diff_backend`.
-
-Only valid for manifolds that are embedded in a special way in the Euclidean space.
-See [^Absil2008], Section 3.6.1 for details.
-
-[^Absil2008]:
-    > Absil, P. A., et al. Optimization Algorithms on Matrix Manifolds. 2008.
-"""
-struct RiemannianProjectionGradientBackend{TADBackend<:AbstractDiffBackend} <:
+@doc raw"""
+    RiemannianProjectionBackend <: AbstractRiemannianDiffBackend
+
+This backend computes the differentiation in the embedding, which is currently limited
+to the gradient. Let ``mathcal M`` denote a manifold embedded in some ``R^m``, where ``m``
+is usually (much) larger than the manifold dimension.
+Then we require three tools
+
+* A function ``f̃: ℝ^m → ℝ`` such that its restriction to the manifold yields the cost
+  function ``f`` of interest.
+* A [`project`](@ref) function to project tangent vectors from the embedding (at `T_pℝ^m``)
+  back onto the tangent space ``T_p\mathcal M``. This also includes possible changes
+  of the representation of the tangent vector (e.g. in the Lie algebra or in a different data format).
+* A [`change_representer`](@ref) for non-isometrically embedded manifolds,
+  i.e. where the tangent space ``T_p\mathcal M`` of the manifold does not inherit
+  the inner product from restriction of the inner product from the tangent space ``T_pℝ^m``
+  of the embedding
+
+For more details see [^AbsilMahonySepulchre2008], Section 3.6.1 for a derivation on submanifolds.
+
+[^AbsilMahonySepulchre2008]:
+    > Absil, P.-A., Mahony, R. and Sepulchre R.,
+    > _Optimization Algorithms on Matrix Manifolds_
+    > Princeton University Press, 2008,
+    > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244)
+    > [open access](http://press.princeton.edu/chapters/absil/)
+"""
+struct RiemannianProjectionBackend{TADBackend<:AbstractDiffBackend} <:
        AbstractRiemannianDiffBackend
     diff_backend::TADBackend
 end
 
-function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionGradientBackend)
+function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionBackend)
     amb_grad = _gradient(f, p, backend.diff_backend)
     return change_representer(M, EuclideanMetric(), p, project(M, p, amb_grad))
 end
 
-function gradient!(
-    M::AbstractManifold,
-    f,
-    X,
-    p,
-    backend::RiemannianProjectionGradientBackend,
-)
+function gradient!(M::AbstractManifold, f, X, p, backend::RiemannianProjectionBackend)
     amb_grad = embed(M, p, X)
     _gradient!(f, amb_grad, p, backend.diff_backend)
     project!(M, X, p, amb_grad)
diff --git a/src/differentiation/zygote.jl b/src/differentiation/zygote.jl
index ffe180113f..883e027e30 100644
--- a/src/differentiation/zygote.jl
+++ b/src/differentiation/zygote.jl
@@ -8,4 +8,6 @@ function Manifolds._gradient!(f, X, p, ::ZygoteDiffBackend)
     return copyto!(X, Zygote.gradient(f, p)[1])
 end
 
-push!(Manifolds._diff_backends, ZygoteDiffBackend())
+if default_differential_backend() === NoneDiffBackend()
+    set_default_differential_backend!(ZygoteDiffBackend())
+end
diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl
index 38dd212a28..57d9c39bcc 100644
--- a/src/manifolds/ConnectionManifold.jl
+++ b/src/manifolds/ConnectionManifold.jl
@@ -59,7 +59,7 @@ end
         M::AbstractManifold,
         p,
         B::AbstractBasis;
-        backend::AbstractDiffBackend = diff_backend(),
+        backend::AbstractDiffBackend = default_differential_backend(),
     )
 
 Compute the Christoffel symbols of the second kind in local coordinates of basis `B`.
@@ -74,6 +74,7 @@ where ``Γ_{ijk}`` are the Christoffel symbols of the first kind
 (see [`christoffel_symbols_first`](@ref)), and ``g^{kl}`` is the inverse of the local
 representation of the metric tensor. The dimensions of the resulting multi-dimensional array
 are ordered ``(l,i,j)``.
+
 """
 christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis)
 
@@ -89,7 +90,7 @@ christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis)
         M::AbstractManifold,
         p,
         B::AbstractBasis;
-        backend::AbstractDiffBackend = diff_backend(),
+        backend::AbstractDiffBackend = default_differential_backend(),
     )
 
 Get partial derivatives of the Christoffel symbols of the second kind
@@ -106,15 +107,25 @@ function christoffel_symbols_second_jacobian(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
-    n = size(p, 1)
+    d = manifold_dimension(M)
     ∂Γ = reshape(
-        _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend),
-        n,
-        n,
-        n,
-        n,
+        _jacobian(
+            c -> christoffel_symbols_second(
+                M,
+                retract(M, p, get_vector(M, p, c, B), retraction),
+                B;
+                backend=backend,
+            ),
+            p,
+            backend,
+        ),
+        d,
+        d,
+        d,
+        d,
     )
     return ∂Γ
 end
@@ -152,17 +163,17 @@ in an embedded space.
 exp(::AbstractConnectionManifold, ::Any...)
 
 @decorator_transparent_fallback function exp!(M::AbstractConnectionManifold, q, p, X)
-    tspan = (0.0, 1.0)
-    A = get_default_atlas(M)
-    i = get_chart_index(M, A, p)
-    B = induced_basis(M, A, i, TangentSpace)
-    sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0])
-    n = length(p)
-    return copyto!(q, sol.u[1][(n + 1):end])
+    return retract!(
+        M,
+        q,
+        p,
+        X,
+        ODEExponentialRetraction(ManifoldsBase.default_retraction_method(M)),
+    )
 end
 
 """
-    gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())
+    gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())
 
 Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`.
 This is equal to half of the scalar Ricci curvature, see [`ricci_curvature`](@ref).
@@ -195,7 +206,7 @@ function injectivity_radius(M::AbstractConnectionManifold, p, m::ExponentialRetr
 end
 
 """
-    ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())
+    ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())
 
 Compute the Ricci tensor, also known as the Ricci curvature tensor,
 of the manifold `M` at the point `p` using basis `B`,
@@ -217,7 +228,7 @@ end
 )
 
 @doc raw"""
-    riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=diff_backend())
+    riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend())
 
 Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature
 tensor, at the point `p` in local coordinates defined by `B`. The dimensions of the
@@ -236,11 +247,19 @@ function riemann_tensor(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
     n = size(p, 1)
-    Γ = christoffel_symbols_second(M, p, B; backend=backend)
-    ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n
+    Γ = christoffel_symbols_second(M, p, B; backend=backend, retraction=retraction)
+    ∂Γ =
+        christoffel_symbols_second_jacobian(
+            M,
+            p,
+            B;
+            backend=backend,
+            retraction=retraction,
+        ) ./ n
     R = allocate(∂Γ, Size(n, n, n, n))
     @einsum R[l, i, j, k] =
         ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k]
@@ -255,17 +274,17 @@ end
 
 @doc raw"""
     solve_exp_ode(
-        M::AbstractConnectionManifold,
+        M::AbstractManifold,
         p,
         X,
-        tspan,
         B::AbstractBasis;
-        backend::AbstractDiffBackend = diff_backend(),
+        backend::AbstractDiffBackend = default_differential_backend(),
         solver = AutoVern9(Rodas5()),
+        retraction = ManifoldsBase.default_retraction_method(M),
         kwargs...,
     )
 
-Approximate the exponential map on the manifold over the provided timespan
+Approximate the exponential map on the manifold by evaluating the ODE descripting the geodesic at 1,
 assuming the default connection of the given manifold by solving the ordinary differential
 equation
 
@@ -278,9 +297,8 @@ the Einstein summation convention is assumed. The arguments `tspan` and
 `solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword
 arguments that will be passed to `OrdinaryDiffEq.solve`.
 
-Currently, the numerical integration is only accurate when using a single
-coordinate chart that covers the entire manifold. This excludes coordinates
-in an embedded space.
+Using a retraction and a basis, the solver works on any intrinsic manifold that
+provides a basis for tangent spaces around the start point `p`.
 
 !!! note
     This function only works when
@@ -289,9 +307,9 @@ in an embedded space.
     using OrdinaryDiffEq
     ```
 """
-function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...)
+function solve_exp_ode(M, p, X; kwargs...)
     return error(
-        "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.",
+        "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.",
     )
 end
 
diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl
index 71888e6133..f542ee79fc 100644
--- a/src/manifolds/Euclidean.jl
+++ b/src/manifolds/Euclidean.jl
@@ -199,7 +199,16 @@ function get_basis(M::Euclidean, p, B::DiagonalizingOrthonormalBasis)
     return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, eigenvalues, vecs))
 end
 
-function get_coordinates!(M::Euclidean, Y, p, X, ::DefaultOrDiagonalizingBasis{ℝ})
+function get_coordinates!(
+    M::Euclidean,
+    Y,
+    p,
+    X,
+    ::Union{
+        DefaultOrDiagonalizingBasis{ℝ},
+        InducedBasis{ℝ,TangentSpaceType,<:RetractionAtlas},
+    },
+)
     S = representation_size(M)
     PS = prod(S)
     copyto!(Y, reshape(X, PS))
@@ -218,7 +227,16 @@ function get_coordinates!(
     return Y
 end
 
-function get_vector!(M::Euclidean, Y, ::Any, X, ::DefaultOrDiagonalizingBasis{ℝ})
+function get_vector!(
+    M::Euclidean,
+    Y,
+    ::Any,
+    X,
+    ::Union{
+        DefaultOrDiagonalizingBasis{ℝ},
+        InducedBasis{ℝ,TangentSpaceType,<:RetractionAtlas},
+    },
+)
     S = representation_size(M)
     copyto!(Y, reshape(X, S))
     return Y
@@ -239,7 +257,6 @@ function get_vector!(M::Euclidean{<:Tuple,ℂ}, Y, ::Any, X, ::DefaultOrDiagonal
     copyto!(Y, reshape(X[1:N] + im * X[(N + 1):end], S))
     return Y
 end
-
 @doc raw"""
     injectivity_radius(M::Euclidean)
 
diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl
index 2eb7db887b..95bacabafd 100644
--- a/src/manifolds/MetricManifold.jl
+++ b/src/manifolds/MetricManifold.jl
@@ -200,7 +200,7 @@ end
         M::MetricManifold,
         p,
         B::AbstractBasis;
-        backend::AbstractDiffBackend = diff_backend(),
+        backend::AbstractDiffBackend = default_differential_backend(),
     )
 
 Compute the Christoffel symbols of the first kind in local coordinates of basis `B`.
@@ -219,9 +219,10 @@ function christoffel_symbols_first(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
-    ∂g = local_metric_jacobian(M, p, B; backend=backend)
+    ∂g = local_metric_jacobian(M, p, B; backend=backend, retraction=retraction)
     n = size(∂g, 1)
     Γ = allocate(∂g, Size(n, n, n))
     @einsum Γ[i, j, k] = 1 / 2 * (∂g[k, j, i] + ∂g[i, k, j] - ∂g[i, j, k])
@@ -238,10 +239,11 @@ function christoffel_symbols_second(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
     Ginv = inverse_local_metric(M, p, B)
-    Γ₁ = christoffel_symbols_first(M, p, B; backend=backend)
+    Γ₁ = christoffel_symbols_first(M, p, B; backend=backend, retraction=retraction)
     Γ₂ = allocate(Γ₁)
     @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k]
     return Γ₂
@@ -272,7 +274,7 @@ end
     B::AbstractBasis,
 )
 """
-    einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())
+    einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_badefault_differential_backendckend())
 
 Compute the Einstein tensor of the manifold `M` at the point `p`, see [https://en.wikipedia.org/wiki/Einstein_tensor](https://en.wikipedia.org/wiki/Einstein_tensor)
 """
@@ -281,7 +283,7 @@ function einstein_tensor(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
 )
     Ric = ricci_tensor(M, p, B; backend=backend)
     g = local_metric(M, p, B)
@@ -452,22 +454,32 @@ local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
         M::AbstractManifold,
         p,
         B::AbstractBasis;
-        backend::AbstractDiffBackend = diff_backend(),
+        backend::AbstractDiffBackend,
     )
 
 Get partial derivatives of the local metric of `M` at `p` in basis `B` with respect to the
 coordinates of `p`, ``\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}``. The
 dimensions of the resulting multi-dimensional array are ordered ``(i,j,k)``.
 """
-local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
+local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis, ::AbstractDiffBackend)
 function local_metric_jacobian(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
-    n = size(p, 1)
-    ∂g = reshape(_jacobian(q -> local_metric(M, q, B), p, backend), n, n, n)
+    d = manifold_dimension(M)
+    ∂g = reshape(
+        _jacobian(
+            c -> local_metric(M, retract(M, p, get_vector(M, p, c, B), retraction), B),
+            zeros(d),
+            backend,
+        ),
+        d,
+        d,
+        d,
+    )
     return ∂g
 end
 @decorator_transparent_signature local_metric_jacobian(
@@ -516,7 +528,7 @@ function metric(M::MetricManifold)
 end
 
 @doc raw"""
-    ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())
+    ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())
 
 Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`.
 The curvature is computed as the trace of the Ricci curvature tensor with respect to
@@ -531,10 +543,11 @@ function ricci_curvature(
     M::AbstractManifold,
     p,
     B::AbstractBasis;
-    backend::AbstractDiffBackend=diff_backend(),
+    backend::AbstractDiffBackend=default_differential_backend(),
+    retraction::AbstractRetractionMethod=ManifoldsBase.default_retraction_method(M),
 )
     Ginv = inverse_local_metric(M, p, B)
-    Ric = ricci_tensor(M, p, B; backend=backend)
+    Ric = ricci_tensor(M, p, B; backend=backend, retraction=retraction)
     S = sum(Ginv .* Ric)
     return S
 end
diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl
index 92f07d1a25..669e5886a0 100644
--- a/src/manifolds/Sphere.jl
+++ b/src/manifolds/Sphere.jl
@@ -308,6 +308,17 @@ function inverse_retract!(::AbstractSphere, X, p, q, ::ProjectionInverseRetracti
     return (X .= q ./ real(dot(p, q)) .- p)
 end
 
+@doc raw"""
+    local_metric(M::Sphere{n}, p, ::DefaultOrthonormalBasis)
+
+return the local representation of the metric in a [`DefaultOrthonormalBasis`](@ref), namely
+the diagonal matrix of size ``n×n`` with ones on the diagonal, since the metric is obtained
+from the embedding by restriction to the tangent space ``T_p\mathcal M`` at ``p``.
+"""
+function local_metric(::Sphere{n,ℝ}, p, B::DefaultOrthonormalBasis) where {n}
+    return Diagonal(ones(SVector{n,eltype(p)}))
+end
+
 @doc raw"""
     log(M::AbstractSphere, p, q)
 
diff --git a/src/ode.jl b/src/ode.jl
deleted file mode 100644
index 48433b05fc..0000000000
--- a/src/ode.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-function solve_exp_ode(
-    M::MetricManifold,
-    x,
-    v,
-    tspan,
-    B::AbstractBasis;
-    solver=AutoVern9(Rodas5()),
-    backend=diff_backend(),
-    kwargs...,
-)
-    n = length(x)
-    iv = SVector{n}(1:n)
-    ix = SVector{n}((n + 1):(2n))
-    u0 = allocate(x, 2n)
-    u0[iv] .= v
-    u0[ix] .= x
-
-    function exp_problem(u, p, t)
-        M = p[1]
-        dx = u[iv]
-        x = u[ix]
-        ddx = allocate(u, Size(n))
-        du = allocate(u)
-        Γ = christoffel_symbols_second(M, x, B; backend=backend)
-        @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j]
-        du[iv] .= ddx
-        du[ix] .= dx
-        return Base.convert(typeof(u), du)
-    end
-
-    p = (M,)
-    prob = ODEProblem(exp_problem, u0, tspan, p)
-    sol = solve(prob, solver; kwargs...)
-    return sol
-end
diff --git a/test/ambiguities.jl b/test/ambiguities.jl
index 57988178d2..6377f04ccf 100644
--- a/test/ambiguities.jl
+++ b/test/ambiguities.jl
@@ -4,7 +4,7 @@
         # Interims solution until we follow what was proposed in
         # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
         fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs)
-        FMBS_LIMIT = 20
+        FMBS_LIMIT = 21
         @test length(fmbs) <= FMBS_LIMIT
         if length(fmbs) > FMBS_LIMIT
             for amb in fmbs
@@ -16,7 +16,7 @@
         # Interims solution until we follow what was proposed in
         # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
         fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms)
-        FMS_LIMIT = 21
+        FMS_LIMIT = 22
         if length(fms) > FMS_LIMIT
             for amb in fms
                 println(amb)
diff --git a/test/differentiation.jl b/test/differentiation.jl
index cc0d2361ed..9366dd2f00 100644
--- a/test/differentiation.jl
+++ b/test/differentiation.jl
@@ -10,99 +10,76 @@ using Manifolds:
     _gradient,
     _gradient!
 
-using FiniteDifferences
+using FiniteDifferences, FiniteDiff
 using LinearAlgebra: Diagonal, dot
 
 @testset "Differentiation backend" begin
     fd51 = Manifolds.FiniteDifferencesBackend()
-    @testset "diff_backend" begin
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-        @test length(diff_backends()) == 3
-        @test diff_backends()[1] isa Manifolds.FiniteDifferencesBackend
+    @testset "default_differential_backend" begin
+        #ForwardDiff is loaded first in utils.
+        @test default_differential_backend() === Manifolds.ForwardDiffBackend()
 
         @test length(fd51.method.grid) == 5
         # check method order
         @test typeof(fd51.method).parameters[2] == 1
         fd71 = Manifolds.FiniteDifferencesBackend(central_fdm(7, 1))
-        @test diff_backend!(fd71) == fd71
-        @test diff_backend() == fd71
+        @test set_default_differential_backend!(fd71) == fd71
+        @test default_differential_backend() == fd71
     end
 
     using ForwardDiff
 
     fwd_diff = Manifolds.ForwardDiffBackend()
     @testset "ForwardDiff" begin
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-        @test length(diff_backends()) == 3
-        @test diff_backends()[1] isa Manifolds.FiniteDifferencesBackend
-        @test diff_backends()[2] == fwd_diff
-
-        @test diff_backend!(fwd_diff) == fwd_diff
-        @test diff_backend() == fwd_diff
-        @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-
-        diff_backend!(fwd_diff)
-        @test diff_backend() == fwd_diff
-        diff_backend!(fd51)
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
+
+        @test set_default_differential_backend!(fwd_diff) == fwd_diff
+        @test default_differential_backend() == fwd_diff
+        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
+
+        set_default_differential_backend!(fwd_diff)
+        @test default_differential_backend() == fwd_diff
+        set_default_differential_backend!(fd51)
     end
 
     using FiniteDiff
 
     finite_diff = Manifolds.FiniteDiffBackend()
     @testset "FiniteDiff" begin
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-        @test length(diff_backends()) == 4
-        @test diff_backends()[4] == finite_diff
-
-        @test diff_backend!(finite_diff) == finite_diff
-        @test diff_backend() == finite_diff
-        @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-
-        diff_backend!(finite_diff)
-        @test diff_backend() == finite_diff
-        diff_backend!(fd51)
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
+
+        @test set_default_differential_backend!(finite_diff) == finite_diff
+        @test default_differential_backend() == finite_diff
+        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
+
+        set_default_differential_backend!(finite_diff)
+        @test default_differential_backend() == finite_diff
+        set_default_differential_backend!(fd51)
     end
 
     using ReverseDiff
 
     reverse_diff = Manifolds.ReverseDiffBackend()
     @testset "ReverseDiff" begin
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-        @test length(diff_backends()) == 4
-        @test diff_backends()[3] == reverse_diff
-
-        @test diff_backend!(reverse_diff) == reverse_diff
-        @test diff_backend() == reverse_diff
-        @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-
-        diff_backend!(reverse_diff)
-        @test diff_backend() == reverse_diff
-        diff_backend!(fd51)
-    end
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
 
-    using Zygote: Zygote
+        @test set_default_differential_backend!(reverse_diff) == reverse_diff
+        @test default_differential_backend() == reverse_diff
+        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
+        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend
 
-    zygote_diff = Manifolds.ZygoteDiffBackend()
-    @testset "Zygote" begin
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-        @test length(diff_backends()) == 5
-        @test diff_backends()[5] == zygote_diff
-
-        @test diff_backend!(zygote_diff) == zygote_diff
-        @test diff_backend() == zygote_diff
-        @test diff_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
-        @test diff_backend() isa Manifolds.FiniteDifferencesBackend
-
-        diff_backend!(zygote_diff)
-        @test diff_backend() == zygote_diff
-        diff_backend!(fd51)
+        set_default_differential_backend!(reverse_diff)
+        @test default_differential_backend() == reverse_diff
+        set_default_differential_backend!(fd51)
     end
 
+    using Zygote
+    zygote_diff = Manifolds.ZygoteDiffBackend()
+
     @testset "gradient" begin
-        diff_backend!(fd51)
+        set_default_differential_backend!(fd51)
         r2 = Euclidean(2)
 
         c1(t) = [sin(t), cos(t)]
@@ -121,59 +98,59 @@ using LinearAlgebra: Diagonal, dot
         end
 
         @testset for backend in [fd51, fwd_diff, finite_diff]
-            diff_backend!(backend)
+            set_default_differential_backend!(backend)
             @test _derivative(c1, 0.0) ≈ [1.0, 0.0]
             X = [-1.0, -1.0]
             @test _derivative!(c1, X, 0.0) === X
             @test isapprox(X, [1.0, 0.0])
         end
         @testset for backend in [fd51, fwd_diff, finite_diff, reverse_diff, zygote_diff]
-            diff_backend!(backend)
+            set_default_differential_backend!(backend)
             X = [-1.0, -1.0]
             @test _gradient(f1, [1.0, -1.0]) ≈ [1.0, -2.0]
             @test _gradient!(f1, X, [1.0, -1.0]) === X
             @test X ≈ [1.0, -2.0]
         end
-        diff_backend!(Manifolds.NoneDiffBackend())
+        set_default_differential_backend!(Manifolds.NoneDiffBackend())
         @testset for backend in [fd51, Manifolds.ForwardDiffBackend()]
             @test _derivative(c1, 0.0, backend) ≈ [1.0, 0.0]
             @test _gradient(f1, [1.0, -1.0], backend) ≈ [1.0, -2.0]
         end
 
-        diff_backend!(fd51)
+        set_default_differential_backend!(fd51)
     end
 end
 
-rb_onb_default = RiemannianONBDiffBackend(
-    diff_backend(),
+rb_onb_default = TangentDiffBackend(
+    default_differential_backend(),
     Manifolds.ExponentialRetraction(),
     Manifolds.LogarithmicInverseRetraction(),
     DefaultOrthonormalBasis(),
 )
 
-rb_onb_fd51 = RiemannianONBDiffBackend(
+rb_onb_fd51 = TangentDiffBackend(
     Manifolds.FiniteDifferencesBackend(),
     Manifolds.ExponentialRetraction(),
     Manifolds.LogarithmicInverseRetraction(),
     DefaultOrthonormalBasis(),
 )
 
-rb_onb_fwd_diff = RiemannianONBDiffBackend(
+rb_onb_fwd_diff = TangentDiffBackend(
     Manifolds.ForwardDiffBackend(),
     Manifolds.ExponentialRetraction(),
     Manifolds.LogarithmicInverseRetraction(),
     DefaultOrthonormalBasis(),
 )
 
-rb_onb_finite_diff = RiemannianONBDiffBackend(
+rb_onb_finite_diff = TangentDiffBackend(
     Manifolds.FiniteDiffBackend(),
     Manifolds.ExponentialRetraction(),
     Manifolds.LogarithmicInverseRetraction(),
     DefaultOrthonormalBasis(),
 )
 
-rb_onb_default2 = RiemannianONBDiffBackend(
-    diff_backend(),
+rb_onb_default2 = TangentDiffBackend(
+    default_differential_backend(),
     Manifolds.ExponentialRetraction(),
     Manifolds.LogarithmicInverseRetraction(),
     CachedBasis(
@@ -182,23 +159,7 @@ rb_onb_default2 = RiemannianONBDiffBackend(
     ),
 )
 
-rb_proj = Manifolds.RiemannianProjectionGradientBackend(diff_backend())
-
-@testset "rdiff_ functions" begin
-    @test Manifolds.rdifferential_backend() ===
-          Manifolds._current_rdifferential_backend.backend
-    @test Manifolds.rgradient_backend() === Manifolds._current_rgradient_backend.backend
-
-    tmp_diff = Manifolds.rdifferential_backend()
-    Manifolds.rdifferential_backend!(rb_onb_finite_diff)
-    @test Manifolds.rdifferential_backend() === rb_onb_finite_diff
-    Manifolds.rdifferential_backend!(tmp_diff)
-
-    tmp_grad = Manifolds.rgradient_backend()
-    Manifolds.rgradient_backend!(rb_onb_finite_diff)
-    @test Manifolds.rgradient_backend() === rb_onb_finite_diff
-    Manifolds.rgradient_backend!(tmp_grad)
-end
+rb_proj = Manifolds.RiemannianProjectionBackend(default_differential_backend())
 
 @testset "Riemannian differentials" begin
     s2 = Sphere(2)
@@ -207,9 +168,9 @@ end
     c1(t) = geodesic(s2, q, p, t)
 
     Xval = [-sqrt(2) / 2, 0.0, sqrt(2) / 2]
-    @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4), Xval)
+    @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4, rb_onb_default), Xval)
     X = similar(p)
-    differential!(s2, c1, X, π / 4)
+    differential!(s2, c1, X, π / 4, rb_onb_default)
     @test isapprox(s2, c1(π / 4), X, Xval)
 
     @testset for backend in [rb_onb_fd51, rb_onb_fwd_diff, rb_onb_finite_diff]
@@ -225,13 +186,10 @@ end
     f1(p) = p[1]
 
     q = [sqrt(2) / 2, 0, sqrt(2) / 2]
-    @test isapprox(s2, q, gradient(s2, f1, q), [0.5, 0.0, -0.5])
     for backend in [rb_onb_default, rb_proj]
         @test isapprox(s2, q, gradient(s2, f1, q, backend), [0.5, 0.0, -0.5])
     end
     X = similar(q)
-    gradient!(s2, f1, X, q)
-    @test isapprox(s2, q, X, [0.5, 0.0, -0.5])
     for backend in [rb_onb_default, rb_proj]
         gradient!(s2, f1, X, q, backend)
         @test isapprox(s2, q, X, [0.5, 0.0, -0.5])
diff --git a/test/metric.jl b/test/metric.jl
index a1562bfe2c..a08101da42 100644
--- a/test/metric.jl
+++ b/test/metric.jl
@@ -1,13 +1,23 @@
 using FiniteDifferences, ForwardDiff
 using LinearAlgebra: I
 using StatsBase: AbstractWeights, pweights
-import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection
+import Manifolds:
+    mean!, median!, InducedBasis, induced_basis, get_chart_index, connection, retract!
+import ManifoldsBase: default_retraction_method
 
 include("utils.jl")
 
 struct TestEuclidean{N} <: AbstractManifold{ℝ} end
 struct TestEuclideanMetric <: AbstractMetric end
 struct TestScaledEuclideanMetric <: AbstractMetric end
+struct TestRetraction <: AbstractRetractionMethod end
+
+ManifoldsBase.default_retraction_method(::TestEuclidean) = TestRetraction()
+function ManifoldsBase.default_retraction_method(
+    ::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric},
+)
+    return TestRetraction()
+end
 
 Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N
 function Manifolds.local_metric(
@@ -36,7 +46,10 @@ function Manifolds.get_coordinates!(
     c,
     ::Any,
     X,
-    ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+    ::Union{
+        DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+        DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+    },
 )
     c .= 1 ./ [1.0:manifold_dimension(M)...] .* X
     return c
@@ -46,7 +59,11 @@ function Manifolds.get_vector!(
     X,
     ::Any,
     c,
-    ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+    ::Union{
+        DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+        DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+        InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+    },
 )
     X .= [1.0:manifold_dimension(M)...] .* c
     return X
@@ -56,7 +73,7 @@ function Manifolds.get_coordinates!(
     c,
     ::Any,
     X,
-    ::DefaultOrthogonalBasis,
+    ::DefaultOrthonormalBasis,
 )
     c .= 1 ./ (2 .* [1.0:manifold_dimension(M)...]) .* X
     return c
@@ -66,22 +83,49 @@ function Manifolds.get_vector!(
     X,
     ::Any,
     c,
-    ::DefaultOrthogonalBasis,
+    ::DefaultOrthonormalBasis,
 )
     X .= 2 .* [1.0:manifold_dimension(M)...] .* c
     return X
 end
+function retract!(
+    ::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric},
+    q,
+    p,
+    X,
+    ::TestRetraction,
+)
+    copyto!(q, p + X)
+    return q
+end
+
 struct TestSphere{N,T} <: AbstractManifold{ℝ}
     r::T
 end
 
 struct TestSphericalMetric <: AbstractMetric end
 
+function retract!(
+    ::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric},
+    q,
+    p,
+    X,
+    ::ProjectionRetraction,
+)
+    copyto!(q, (p + X) ./ norm(p + X))
+    return q
+end
+function ManifoldsBase.default_retraction_method(
+    ::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric},
+)
+    return ProjectionRetraction()
+end
+
 Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N
 function Manifolds.local_metric(
     M::MetricManifold{ℝ,<:TestSphere,<:TestSphericalMetric},
     p,
-    ::InducedBasis,
+    ::Union{InducedBasis,DefaultOrthonormalBasis},
 )
     r = base_manifold(M).r
     d = allocate(p)
@@ -89,6 +133,19 @@ function Manifolds.local_metric(
     d[2] = d[1] * sin(p[1])^2
     return Diagonal(d)
 end
+function Manifolds.get_vector!(
+    ::MetricManifold{ℝ,<:TestSphere{N},<:TestSphericalMetric},
+    Y,
+    p,
+    c,
+    ::Union{
+        DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+        InducedBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+    },
+) where {N}
+    Y .= 1
+    return Y # this is just a dummy to check that dispatch works
+end
 sph_to_cart(θ, ϕ) = [cos(ϕ) * sin(θ), sin(ϕ) * sin(θ), cos(θ)]
 
 struct BaseManifold{N} <: AbstractManifold{ℝ} end
@@ -216,17 +273,6 @@ end
         @test length(methods(is_default_metric)) == 2
     end
 
-    @testset "solve_exp_ode error message" begin
-        E = TestEuclidean{3}()
-        g = TestEuclideanMetric()
-        M = MetricManifold(E, g)
-
-        p = [1.0, 2.0, 3.0]
-        X = [2.0, 3.0, 4.0]
-        @test_throws ErrorException exp(M, p, X)
-        using OrdinaryDiffEq
-        exp(M, p, X)
-    end
     @testset "Local Metric Error message" begin
         M = MetricManifold(BaseManifold{2}(), NotImplementedMetric())
         A = Manifolds.get_default_atlas(M)
@@ -256,8 +302,13 @@ end
 
         i_zeros = get_chart_index(M, A, zeros(3))
         B_i_zeros = induced_basis(M, A, i_zeros, TangentSpace)
-        @test_throws MethodError local_metric_jacobian(E, zeros(3), B_i_zeros)
-        @test_throws MethodError christoffel_symbols_second_jacobian(E, zeros(3), B_i_zeros)
+        # get_vector! not implemented -> ErrorException
+        @test_throws ErrorException local_metric_jacobian(E, zeros(3), B_i_zeros)
+        @test_throws ErrorException christoffel_symbols_second_jacobian(
+            E,
+            zeros(3),
+            B_i_zeros,
+        )
 
         for vtype in (Vector, MVector{n})
             p, X, Y = vtype(randn(n)), vtype(randn(n)), vtype(randn(n))
@@ -278,11 +329,6 @@ end
             @test inner(M, p, fX, fY) ≈ dot(X, G * Y) atol = 1e-6
             @test norm(M, p, fX) ≈ sqrt(dot(X, G * X)) atol = 1e-6
 
-            if VERSION ≥ v"1.1"
-                T = 0:0.5:10
-                @test geodesic(M, p, X, T) ≈ [p + t * X for t in T] atol = 1e-6
-            end
-
             @test christoffel_symbols_first(M, p, B_chart_p) ≈ zeros(n, n, n) atol = 1e-6
             @test christoffel_symbols_second(M, p, B_chart_p) ≈ zeros(n, n, n) atol = 1e-6
             @test riemann_tensor(M, p, B_chart_p) ≈ zeros(n, n, n, n) atol = 1e-6
@@ -356,18 +402,6 @@ end
                 -sin(θ) 0
             ] * X
 
-            if !Sys.iswindows() || Sys.ARCH == :x86_64
-                @testset "numerically integrated geodesics for $vtype" begin
-                    T = 0:0.1:1
-                    @test isapprox(
-                        [sph_to_cart(yi...) for yi in geodesic(M, p, X, T)],
-                        geodesic(S, pcart, Xcart, T);
-                        atol=1e-3,
-                        rtol=1e-3,
-                    )
-                end
-            end
-
             Γ₁ = christoffel_symbols_first(M, p, B_p)
             for i in 1:n, j in 1:n, k in 1:n
                 if (i, j, k) == (1, 2, 2) || (i, j, k) == (2, 1, 2)
@@ -500,8 +534,9 @@ end
         chart_p = get_chart_index(MM2, A, p)
         B_p = induced_basis(MM2, A, chart_p, TangentSpace)
         @test_throws MethodError local_metric(MM2, p, B_p)
-        @test_throws MethodError local_metric_jacobian(MM2, p, B_p)
-        @test_throws MethodError christoffel_symbols_second_jacobian(MM2, p, B_p)
+        # the following two default to solving the ODE, but fail since get_vector is not implemented
+        @test_throws ErrorException local_metric_jacobian(MM2, p, B_p)
+        @test_throws ErrorException christoffel_symbols_second_jacobian(MM2, p, B_p)
         # MM falls back to nondefault error
         @test_throws MethodError projected_distribution(MM, 1, p)
         @test_throws MethodError projected_distribution(MM, 1)
diff --git a/test/ode.jl b/test/ode.jl
new file mode 100644
index 0000000000..49897336a8
--- /dev/null
+++ b/test/ode.jl
@@ -0,0 +1,69 @@
+include("utils.jl")
+
+using FiniteDifferences, ForwardDiff
+using LinearAlgebra: I
+import Manifolds: retract!
+import ManifoldsBase: manifold_dimension, default_retraction_method
+
+#
+# Part I: Euclidean
+#
+struct TestODEEuclidean{N} <: AbstractManifold{ℝ} end
+struct TestODEEuclideanMetric <: AbstractMetric end
+#
+# Part II Spherical
+#
+struct TestODESphere{N} <: AbstractManifold{ℝ} end
+struct TestODESphericalMetric <: AbstractMetric end
+
+function default_retraction_method(
+    ::MetricManifold{ℝ,<:TestODESphere,<:TestODESphericalMetric},
+)
+    return ProjectionRetraction()
+end
+
+manifold_dimension(::TestODESphere{N}) where {N} = N
+function Manifolds.retract!(
+    ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric},
+    q,
+    p,
+    X,
+    ::ProjectionRetraction,
+) where {N}
+    return retract!(Sphere(N), q, p, X)
+end
+function Manifolds.local_metric(
+    ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric},
+    p,
+    B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+) where {N}
+    return Manifolds.local_metric(MetricManifold(Sphere(N), EuclideanMetric()), p, B)
+end
+function Manifolds.get_coordinates!(
+    ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric},
+    c,
+    p,
+    X,
+    B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+) where {N}
+    return get_coordinates!(Sphere(N), c, p, X, B)
+end
+function Manifolds.get_vector!(
+    ::MetricManifold{ℝ,<:TestODESphere{N},<:TestODESphericalMetric},
+    Y,
+    p,
+    c,
+    B::DefaultOrthonormalBasis{ℝ,<:ManifoldsBase.TangentSpaceType},
+) where {N}
+    return get_vector!(Sphere(N), Y, p, c, B)
+end
+#@testset "Test ODE setup for computing geodesics" begin
+M = TestODESphere{2}()
+p = [0.0, 0.0, 1.0]
+X = π / (2 * sqrt(2)) .* [0.0, 1.0, 1.0]
+M2 = MetricManifold(M, TestODESphericalMetric())
+#    @test_throws ErrorException exp(M, p, X)
+#    @test_throws ErrorException exp(M2, p, X)
+using OrdinaryDiffEq
+exp(M2, p, X)
+#end
diff --git a/test/runtests.jl b/test/runtests.jl
index 25e841c377..27115de4e4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -147,6 +147,7 @@ include("utils.jl")
     include_test("manifolds/graph.jl")
 
     include_test("metric.jl")
+    include_test("ode.jl")
     include_test("statistics.jl")
     include_test("approx_inverse_retraction.jl")