diff --git a/analysis/mean_squared_displacement.jl b/analysis/mean_squared_displacement.jl
new file mode 100644
index 0000000..586a547
--- /dev/null
+++ b/analysis/mean_squared_displacement.jl
@@ -0,0 +1,161 @@
+using CairoMakie
+using LaTeXStrings: @L_str
+using Dates: Dates
+using Random: Random
+using StaticArrays: SVector
+using JLD2: JLD2
+
+using ReCo: ReCo
+
+# Disable the periodic boundary conditions
+ReCo.update_verlet_lists!(args...) = nothing
+ReCo.gen_cell_list(args...) = nothing
+ReCo.gen_cell_list_box(args...) = nothing
+ReCo.push_to_verlet_list!(args...) = nothing
+ReCo.restrict_coordinate(value::Float64, args...) = value
+ReCo.restrict_coordinates(v::SVector{2,Float64}, args...) = v
+ReCo.restrict_coordinates!(args...) = nothing
+ReCo.minimum_image_coordinate(value::Float64, args...) = value
+ReCo.minimum_image(v::SVector{2,Float64}, args...) = v
+
+function mean_squared_displacement(;
+    n_simulations::Int64, v₀s::AbstractVector{Float64}, T::Float64
+)
+    Random.seed!(42)
+
+    n_v₀s = length(v₀s)
+
+    δt = ReCo.DEFAULT_δt
+    Dₜ = ReCo.DEFAULT_Dₜ
+
+    main_parent_dir = "mean_squared_displacement_$(Dates.now())"
+    mkpath(main_parent_dir)
+
+    sim_dirs = Matrix{String}(undef, (n_simulations, n_v₀s))
+
+    for (v₀_ind, v₀) in enumerate(v₀s)
+        max_possible_displacement = T * v₀ + T / δt * sqrt(2 * Dₜ * δt)
+
+        parent_dir = main_parent_dir * "/$v₀"
+
+        Threads.@threads for sim_ind in 1:n_simulations
+            dir = ReCo.init_sim(;
+                n_particles=1,
+                v₀=v₀,
+                parent_dir=parent_dir,
+                comment="$sim_ind",
+                half_box_len=max_possible_displacement,
+            )
+
+            sim_dirs[sim_ind, v₀_ind] = dir
+
+            ReCo.run_sim(dir; duration=T, seed=rand(1:typemax(Int64)), snapshot_at=0.01)
+        end
+    end
+
+    ts = Float64[]
+
+    bundle_paths = ReCo.sorted_bundle_paths(sim_dirs[1, 1])
+
+    for i in 2:length(bundle_paths) # Skip the first bundle to avoid t = 0
+        bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
+
+        append!(ts, bundle.t)
+    end
+
+    mean_sq_displacements = zeros((length(ts), n_v₀s))
+
+    @simd for v₀_ind in 1:n_v₀s
+        for sim_ind in 1:n_simulations
+            sim_dir = sim_dirs[sim_ind, v₀_ind]
+
+            bundle_paths = ReCo.sorted_bundle_paths(sim_dir)
+            snapshot_ind = 1
+
+            for i in 2:length(bundle_paths)
+                bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
+
+                for c in bundle.c
+                    sq_displacement = ReCo.sq_norm2d(c)
+                    mean_sq_displacements[snapshot_ind, v₀_ind] += sq_displacement
+
+                    snapshot_ind += 1
+                end
+            end
+        end
+    end
+
+    mean_sq_displacements ./= n_simulations
+
+    return (ts, mean_sq_displacements)
+end
+
+function expected_mean_squared_displacement(t::Float64, v₀::Float64)
+    Dₜ = ReCo.DEFAULT_Dₜ
+    particle_radius = ReCo.DEFAULT_PARTICLE_RADIUS
+
+    Dᵣ = 3 * Dₜ / ((2 * particle_radius)^2)
+
+    return (4 * Dₜ + 2 * v₀^2 / Dᵣ) * t + 2 * v₀^2 * (exp(-Dᵣ * t) - 1) / (Dᵣ^2)
+end
+
+function plot_mean_sq_displacement_with_expectation(
+    ts::Vector{Float64},
+    mean_sq_displacements::Matrix{Float64},
+    v₀s::AbstractVector{Float64},
+)
+    CairoMakie.activate!()
+    set_theme!()
+
+    text_width_in_pt = 405
+
+    fig = Figure(;
+        resolution=(text_width_in_pt, 0.55 * text_width_in_pt),
+        fontsize=10,
+        figure_padding=1,
+    )
+
+    ax = Axis(fig[1, 1]; xlabel=L"t", ylabel=L"\mathbf{MSD}", xscale=log10, yscale=log10)
+
+    t_linrange = LinRange(ts[1], ts[end], 1000)
+
+    v₀_scatter_plots = []
+
+    for (v₀_ind, v₀) in enumerate(v₀s)
+        scatter_plot = scatter!(
+            ax, ts, view(mean_sq_displacements, :, v₀_ind); markersize=4
+        )
+
+        push!(v₀_scatter_plots, scatter_plot)
+
+        expected_mean_sq_displacements = expected_mean_squared_displacement.(t_linrange, v₀)
+
+        lines!(ax, t_linrange, expected_mean_sq_displacements)
+    end
+
+    Legend(fig[1, 2], v₀_scatter_plots, [L"v₀ = %$v₀" for v₀ in v₀s])
+
+    colgap!(fig.layout, 5)
+    rowgap!(fig.layout, 5)
+
+    parent_dir = "exports/graphics/"
+    mkpath(parent_dir)
+
+    save("$parent_dir/mean_squared_displacement.pdf", fig; pt_per_unit=1)
+
+    return nothing
+end
+
+function run_analysis()
+    v₀s = SVector(0.0, 20.0, 40.0, 60.0, 80.0)
+
+    ts, mean_sq_displacements = mean_squared_displacement(;
+        n_simulations=3 * Threads.nthreads(), v₀s=v₀s, T=10.0
+    )
+
+    plot_mean_sq_displacement_with_expectation(ts, mean_sq_displacements, v₀s)
+
+    return nothing
+end
+
+# run_analysis()
\ No newline at end of file