diff --git a/Project.toml b/Project.toml
index 3ca187a..daf63c7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,7 +24,7 @@ LinearAlgebra = "1.6"
 PrecompileTools = "1"
 Printf = "1.6"
 Reexport = "1"
-SatelliteToolboxBase = "0.2, 0.3"
+SatelliteToolboxBase = "0.3.2"
 SatelliteToolboxSgp4 = "2.1"
 SatelliteToolboxTle = "1"
 StaticArrays = "1"
diff --git a/src/SatelliteToolboxPropagators.jl b/src/SatelliteToolboxPropagators.jl
index 85ec27e..766f35f 100644
--- a/src/SatelliteToolboxPropagators.jl
+++ b/src/SatelliteToolboxPropagators.jl
@@ -11,6 +11,8 @@ using StaticArrays
 @reexport using SatelliteToolboxBase
 @reexport using SatelliteToolboxSgp4
 
+import Base: copy
+
 ############################################################################################
 #                                           API                                            #
 ############################################################################################
@@ -47,6 +49,8 @@ include("./api/j4osc.jl")
 include("./api/sgp4.jl")
 include("./api/twobody.jl")
 
+include("./copy.jl")
+
 include("./propagators/j2.jl")
 include("./propagators/j2osc.jl")
 include("./propagators/j4.jl")
diff --git a/src/api/Propagators.jl b/src/api/Propagators.jl
index b6bde3e..bca81e7 100644
--- a/src/api/Propagators.jl
+++ b/src/api/Propagators.jl
@@ -8,15 +8,25 @@ module Propagators
 
 using Dates
 using Crayons
+using StaticArrays
 
-import Base: eltype, length, iterate, show
+import Base: copy, eltype, length, iterate, show
+import SatelliteToolboxBase: @maybe_threads, get_partition
 
 export OrbitPropagator
 
+############################################################################################
+#                                         Contants                                         #
+############################################################################################
+
 # Escape sequences related to the crayons.
 const _D = string(Crayon(reset = true))
 const _B = string(crayon"bold")
 
+############################################################################################
+#                                          Types                                           #
+############################################################################################
+
 """
     abstract type OrbitPropagator{Tepoch<:Number, T<:Number}
 
@@ -24,6 +34,10 @@ Abstract type for the orbit propagators.
 """
 abstract type OrbitPropagator{Tepoch<:Number, T<:Number} end
 
+############################################################################################
+#                                     Public Functions                                     #
+############################################################################################
+
 """
     fit_mean_elements(::Val{:propagator}, vjd::AbstractVector{Tjd}, vr_i::AbstractVector{Tv}, vv_i::AbstractVector{Tv}; kwargs...) where {Tjd<:Number, Tv<:AbstractVector} -> <Mean elements>
 
@@ -127,6 +141,74 @@ Propagate the orbit using `orbp` by `t` [s] from the initial orbit epoch.
 """
 function propagate! end
 
+"""
+    propagate!(orbp::OrbitPropagator{Tepoch, T}, vt::AbstractVector; kwargs...) where {Tepoch <: Number, T <: Number} -> Vector{SVector{3, T}}, Vector{SVector{3, T}}
+
+Propagate the orbit using `orbp` for every instant defined in `vt` [s].
+
+# Keywords
+
+- `ntasks::Integer`: Number of parallel tasks to propagate the orbit. If it is set to a
+    number equal or lower than 1, the function will propagate the orbit sequentially.
+    (**Default** = `Threads.nthreads()`)
+
+# Returns
+
+- `Vector{SVector{3, T}}`: Array with the position vectors [m] in the inertial frame at each
+    propagation instant defined in `vt`.
+- `Vector{SVector{3, T}}`: Array with the velocity vectors [m / s] in the inertial frame at
+    each propagation instant defined in `vt`.
+"""
+function propagate!(
+    orbp::OrbitPropagator{Tepoch, T},
+    vt::AbstractVector;
+    ntasks::Integer = Threads.nthreads()
+) where {Tepoch<:Number, T<:Number}
+    # We need to perform the first propagation to obtain the return type of the propagator.
+    r₀, v₀ = Propagators.propagate!(orbp, first(vt))
+
+    # Number of propagation points.
+    len_vt = length(vt)
+
+    # Allocate the output vectors.
+    vr = Vector{typeof(r₀)}(undef, len_vt)
+    vv = Vector{typeof(v₀)}(undef, len_vt)
+
+    vr[begin] = r₀
+    vv[begin] = v₀
+
+    len_vt == 1 && return vr, vv
+
+    inds = eachindex(vt)
+
+    # We need to store the first index offset of `vt` to allow filling the output vectors
+    # correctly.
+    Δi = firstindex(vt) - 1
+
+    # Make sure the number of tasks is not higher than the number of propagation points.
+    ntasks = min(ntasks, len_vt)
+
+    @maybe_threads ntasks for c in 1:ntasks
+        # We already propagated for the first instant, and we must ensure we propagate the
+        # last instant at the end of the function.
+        i₀, i₁ = get_partition(c, inds[(1 + begin):(end - 1)], ntasks)
+
+        # The propagation usually modifies the structure. Hence we need to copy it for each
+        # task.
+        corbp = c == 1 ? orbp : copy(orbp)
+
+        @inbounds for i in i₀:i₁
+            vr[i - Δi], vv[i - Δi] = Propagators.propagate!(corbp, vt[i])
+        end
+    end
+
+    # We must ensure that the last propagation instant is the obtained at the end to keep
+    # the internal data of the propagation consistent.
+    vr[end], vv[end] = Propagators.propagate!(orbp, last(vt))
+
+    return vr, vv
+end
+
 """
     propagate_to_epoch(::Val{:propagator}, jd::Number, args...; kwargs...)
 
@@ -164,6 +246,28 @@ function propagate_to_epoch!(orbp::OrbitPropagator, jd::Number)
     return propagate!(orbp, 86400 * (jd - epoch(orbp)))
 end
 
+"""
+    propagate_to_epoch!(orbp::OrbitPropagator{Tepoch, T}, vjd::AbstractVector; kwargs...) where {Tepoch, T} -> SVector{3, T}, SVector{3, T}
+
+Propagate the orbit using `orbp` for every epoch defined in `jd` [Julian Day].
+
+# Keywords
+
+- `ntasks::Integer`: Number of parallel tasks to propagate the orbit. If it is set to a
+    number equal or lower than 1, the function will propagate the orbit sequentially.
+    (**Default** = `Threads.nthreads()`)
+
+# Returns
+
+- `Vector{SVector{3, T}}`: Array with the position vectors [m] in the inertial frame at each
+    propagation instant defined in `vt`.
+- `Vector{SVector{3, T}}`: Array with the velocity vectors [m / s] in the inertial frame at
+    each propagation instant defined in `vt`.
+"""
+function propagate_to_epoch!(orbp::OrbitPropagator, vjd::AbstractVector)
+    return propagate!(orbp, 86400 .* (vjd .- epoch(orbp)))
+end
+
 """
     step!(orbp::OrbitPropagator{Tepoch, T}, Δt::Number) where {Tepoch, T} -> SVector{3, T}, SVector{3, T}
 
@@ -180,6 +284,13 @@ function step!(orbp::OrbitPropagator, Δt::Number)
     return propagate!(orbp, last_instant(orbp) + Δt)
 end
 
+############################################################################################
+#                                           Copy                                           #
+############################################################################################
+
+# Define a fallback copy method if the propagator has not implemented one.
+Base.copy(orbp::OrbitPropagator) = deepcopy(orbp)
+
 ############################################################################################
 #                                    Iterator Interface                                    #
 ############################################################################################
diff --git a/src/copy.jl b/src/copy.jl
new file mode 100644
index 0000000..841b600
--- /dev/null
+++ b/src/copy.jl
@@ -0,0 +1,9 @@
+## Description #############################################################################
+#
+# Define the function `copy` for the structures created here.
+#
+############################################################################################
+
+function Base.copy(orbp::OrbitPropagatorSgp4{Tepoch, T}) where {Tepoch <: Number, T <: Number} 
+    return OrbitPropagatorSgp4(copy(orbp.sgp4d))
+end
diff --git a/test/api.jl b/test/api.jl
index 0cf9d4a..5cbdda1 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -67,6 +67,60 @@ struct DummyPropagator{Tepoch, T} <: OrbitPropagator{Tepoch, T} end
     end
 end
 
+@testset "Multi-thread Propagation" verbose = true begin
+    for (T, j2c) in ((Float64, j2c_egm2008), (Float32, j2c_egm2008_f32))
+        @testset "$T" begin
+            jd₀ = date_to_jd(2023, 1, 1, 0, 0, 0)
+            jd₁ = date_to_jd(2023, 1, 5, 0, 0, 0)
+
+            orb = KeplerianElements(
+                jd₀,
+                T(8000e3),
+                T(0.015),
+                T(28.5) |> deg2rad,
+                T(100)  |> deg2rad,
+                T(200)  |> deg2rad,
+                T(45)   |> deg2rad
+            )
+
+            orbp = Propagators.init(Val(:J2), orb; j2c = j2c)
+
+            # == propagate! ================================================================
+
+            ret  = Propagators.propagate!.(orbp, 1:1:100)
+            r, v = Propagators.propagate!(orbp, 1:1:100)
+
+            @test length(r) == 100
+            @test length(v) == 100
+            @test r isa Vector{SVector{3, T}}
+            @test v isa Vector{SVector{3, T}}
+            @test Propagators.last_instant(orbp) == 100.0
+
+            for k in 1:100
+                @test ret[k][1] == r[k]
+                @test ret[k][2] == v[k]
+            end
+
+            # == propagate_to_epoch! =======================================================
+
+            vjd  = collect(jd₀:0.1:jd₁)
+            ret  = Propagators.propagate_to_epoch!.(orbp, vjd)
+            r, v = Propagators.propagate_to_epoch!(orbp, vjd)
+
+            @test length(r) == 41
+            @test length(v) == 41
+            @test r isa Vector{SVector{3, T}}
+            @test v isa Vector{SVector{3, T}}
+            @test Propagators.last_instant(orbp) == 345600.0
+
+            for k in 1:41
+                @test ret[k][1] == r[k]
+                @test ret[k][2] == v[k]
+            end
+        end
+    end
+end
+
 @testset "Show" verbose = true begin
     T = Float64