diff --git a/src/QuanticsTCI.jl b/src/QuanticsTCI.jl
index 0edc347..d1f9ae3 100644
--- a/src/QuanticsTCI.jl
+++ b/src/QuanticsTCI.jl
@@ -12,5 +12,6 @@ export cachedata, quanticsfouriermpo
 
 include("tciinterface.jl")
 include("fouriertransform.jl")
+include("derivative.jl")
 
 end
diff --git a/src/derivative.jl b/src/derivative.jl
new file mode 100644
index 0000000..5ba48e9
--- /dev/null
+++ b/src/derivative.jl
@@ -0,0 +1,115 @@
+# The reference papers for the follwing code are:
+# Quantized tensor networks for solving the Vlasov-Maxwell equations: https://arxiv.org/pdf/2311.07756
+# A quantum-inspired method for solving the Vlasov-Poisson equations: https://journals.aps.org/pre/abstract/10.1103/PhysRevE.106.035208
+
+function S(plus_minus::Symbol)::Matrix{Float64}
+    if plus_minus == :plus
+        return [0 0; 1 0]
+    elseif plus_minus == :minus
+        return [0 1; 0 0]
+    else
+        error("Invalid argument")
+    end
+end
+
+function I()::Matrix{Float64}
+    return [1 0;
+        0 1]
+end
+
+# The first derivative using centered differences with periodic boundary conditions
+function first_order_central_difference_periodic(R::Int64)::Vector{Array{Float64,4}}
+    c1 = 1 / 2
+    t_left = zeros(1, 2, 2, 3)
+    t_left[1, :, :, 1] = I() * 1 / 2^R
+    t_left[1, :, :, 2] = (S(:plus) + S(:minus)) * 1 / 2^R
+    t_left[1, :, :, 3] = (S(:plus) + S(:minus)) * 1 / 2^R
+
+    t_right = zeros(3, 2, 2, 1)
+    t_right[1, :, :, 1] = -c1 * (S(:plus) - S(:minus))
+    t_right[2, :, :, 1] = c1 * S(:plus)
+    t_right[3, :, :, 1] = -c1 * S(:minus)
+
+    t_center = zeros(3, 2, 2, 3)
+    t_center[1, :, :, 1] = I()
+    t_center[1, :, :, 2] = S(:minus)
+    t_center[1, :, :, 3] = S(:plus)
+    t_center[2, :, :, 2] = S(:plus)
+    t_center[3, :, :, 3] = S(:minus)
+
+    mpo = map(1:R) do i
+        if i == 1
+            t_left
+        elseif i == R
+            t_right
+        else
+            t_center
+        end
+    end
+    mpo
+end
+
+# The second-order derivative using centered differences with periodic boundary conditions
+function second_order_central_difference_periodic(R::Int64)::Vector{Array{Float64,4}}
+    c1 = 1.0
+    c0 = -2.0
+    t_left = zeros(1, 2, 2, 3)
+    t_left[1, :, :, 1] = I() * 1 / 2^(2R)
+    t_left[1, :, :, 2] = (S(:plus) + S(:minus)) * 1 / 2^(2R)
+    t_left[1, :, :, 3] = (S(:plus) + S(:minus)) * 1 / 2^(2R)
+
+    t_right = zeros(3, 2, 2, 1)
+    t_right[1, :, :, 1] = (c0 * I() + c1 * (S(:plus) + S(:minus)))
+    t_right[2, :, :, 1] = c1 * S(:minus)
+    t_right[3, :, :, 1] = c1 * S(:plus)
+
+    t_center = zeros(3, 2, 2, 3)
+    t_center[1, :, :, 1] = I()
+    t_center[1, :, :, 2] = S(:plus)
+    t_center[1, :, :, 3] = S(:minus)
+    t_center[2, :, :, 2] = S(:minus)
+    t_center[3, :, :, 3] = S(:plus)
+
+    mpo = map(1:R) do i
+        if i == 1
+            t_left
+        elseif i == R
+            t_right
+        else
+            t_center
+        end
+    end
+    mpo
+end
+
+# linear function av + b on a uniform grid,  where a, b are constants
+# Grid on [−d/2, d/2), N=2^R points, v_i = −d/2 + d_i/R. QTT representation: R cores
+function linear_mps(a::Float64, d::Float64, b::Float64, R::Int64)::Vector{Array{Float64,3}}
+    t_first = zeros(1, 2, 2)
+    t_end = zeros(2, 2, 1)
+    nvec = [0; 1]
+    t_first[1, :, 1] = ones(2)
+    t_first[1, :, 2] = 1 / 2 * a * d * nvec + (-a * d / 2 + b) * ones(2)
+    t_end[1, :, 1] = 1 / 2^R * a * d * nvec
+    t_end[2, :, 1] = ones(2)
+
+    function t_center(i)
+        t_center = zeros(2, 2, 2)
+        t_center[1, :, 1] = ones(2)
+        t_center[2, :, 2] = ones(2)
+        t_center[1, :, 2] = 1 / 2^i * a * d * nvec
+        return t_center
+    end
+
+    mps = map(1:R) do i
+        if i == 1
+            t_first
+        elseif i == R
+            t_end
+        else
+            t_center(i)
+        end
+    end
+
+    return mps
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 2e2ffde..f0be7b3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,3 +7,4 @@ include("test_with_jet.jl")
 include("test_tciinterface.jl")
 include("test_fouriertransform.jl")
 include("test_samplescripts.jl")
+include("test_derivative.jl")
\ No newline at end of file
diff --git a/test/test_derivative.jl b/test/test_derivative.jl
new file mode 100644
index 0000000..aa24f7b
--- /dev/null
+++ b/test/test_derivative.jl
@@ -0,0 +1,81 @@
+import TensorCrossInterpolation as TCI
+import QuanticsTCI as QTCI
+
+function generate_reversed_cartesian_indices(dims)
+    # generate CartesianIndices
+    dims = tuple(dims...)
+    indices = CartesianIndices(dims)
+    # reverse each intex
+    reversed_indices = [CartesianIndex(reverse(Tuple(idx))) for idx in indices]
+    return reversed_indices
+end
+
+function qtt_tomat(tto::TCI.TensorTrain{T,4}) where {T}
+    sitedims = TCI.sitedims(tto)
+    localdims1 = [s[1] for s in sitedims]
+    localdims2 = [s[2] for s in sitedims]
+    mat = Matrix{T}(undef, prod(localdims1), prod(localdims2))
+    for (i, inds1) in enumerate(generate_reversed_cartesian_indices(localdims1)[:])
+        for (j, inds2) in enumerate(generate_reversed_cartesian_indices(localdims2)[:])
+            mat[i, j] = TCI.evaluate(tto, collect(zip(Tuple(inds1), Tuple(inds2))))
+        end
+    end
+    return mat
+end
+
+function qtt_tovec(tt::TCI.TensorTrain{T,3}) where {T}
+    sitedims = TCI.sitedims(tt)
+    localdims1 = [s[1] for s in sitedims]
+    return TCI.evaluate.(Ref(tt), generate_reversed_cartesian_indices(localdims1)[:])
+end
+
+@testset "linear_function" begin
+    mps_linear_ = QTCI.linear_mps(1.0, 8.0, 0.0, 3)
+    mps_linear = TCI.TensorTrain(mps_linear_)
+    @test qtt_tovec(mps_linear) == [-4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0]
+
+end
+
+@testset "first_order_central_difference_periodic" begin
+    R = 3
+    mpo_center_ = QTCI.first_order_central_difference_periodic(R)
+    mpo_center = TCI.TensorTrain(mpo_center_)
+    mat_first_order_central_difference_periodic = (1 / (2 * 2^R)) * [
+        0 1 0 0 0 0 0 -1;
+        -1 0 1 0 0 0 0 0;
+        0 -1 0 1 0 0 0 0;
+        0 0 -1 0 1 0 0 0;
+        0 0 0 -1 0 1 0 0;
+        0 0 0 0 -1 0 1 0;
+        0 0 0 0 0 -1 0 1;
+        1 0 0 0 0 0 -1 0
+    ]
+    @test qtt_tomat(mpo_center) == mat_first_order_central_difference_periodic
+
+    mps_linear_ = QTCI.linear_mps(1.0, 8.0, 0.0, 3)
+    mps_linear = TCI.TensorTrain(mps_linear_)
+    res = TCI.contract(mpo_center, mps_linear) # should be mpo * mps since the order is important
+    qtt_tovec(res) == mat_first_order_central_difference_periodic * qtt_tovec(mps_linear)
+end
+
+@testset "second_order_central_difference_periodic" begin
+    R = 3
+    mpo_center_second_ = QTCI.second_order_central_difference_periodic(R)
+    mpo_center_second = TCI.TensorTrain(mpo_center_second_)
+    mat_second_order_central_difference_periodic = (1 / (2^(2R))) * [
+        -2 1 0 0 0 0 0 1;
+        1 -2 1 0 0 0 0 0;
+        0 1 -2 1 0 0 0 0;
+        0 0 1 -2 1 0 0 0;
+        0 0 0 1 -2 1 0 0;
+        0 0 0 0 1 -2 1 0;
+        0 0 0 0 0 1 -2 1;
+        1 0 0 0 0 0 1 -2
+    ]
+    @test qtt_tomat(mpo_center_second) == mat_second_order_central_difference_periodic
+    mps_linear_ = QTCI.linear_mps(1.0, 8.0, 0.0, 3)
+    mps_linear = TCI.TensorTrain(mps_linear_)
+
+    res = TCI.contract(mpo_center_second, mps_linear) # should be mpo * mps since the order is important
+    qtt_tovec(res) == mat_second_order_central_difference_periodic * qtt_tovec(mps_linear)
+end
\ No newline at end of file