diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index bf293d7d..b5173feb 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -268,10 +268,10 @@ function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: Abstr
 end
 
 """
-    QuadraticSpline(u, t; extrapolate = false)
+    CubicSpline(u, t; extrapolate = false)
 
 It is a spline interpolation using piecewise cubic polynomials between each pair of data points. Its first and second derivative is also continuous.
-Extrapolation extends the last cubic polynomial on each side.
+Second derivative on both ends are zero, which are also called "natural" boundary conditions. Extrapolation extends the last cubic polynomial on each side.
 
 ## Arguments
 
@@ -303,9 +303,9 @@ function CubicSpline(u::uType,
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
-    dl = h[2:(n + 1)]
+    dl = vcat(h[2:n], zero(eltype(h)))
     d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)])
-    du = h[2:(n + 1)]
+    du = vcat(zero(eltype(h)), h[3:(n + 1)])
     tA = Tridiagonal(dl, d_tmp, du)
 
     # zero for element type of d, which we don't know yet
@@ -324,9 +324,9 @@ function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractV
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
-    dl = h[2:(n + 1)]
+    dl = vcat(h[2:n], zero(eltype(h)))
     d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)])
-    du = h[2:(n + 1)]
+    du = vcat(zero(eltype(h)), h[3:(n + 1)])
     tA = Tridiagonal(dl, d_tmp, du)
     d_ = map(
         i -> i == 1 || i == n + 1 ? zeros(eltype(t), size(u[1])) :
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index ebfe1947..5be8c566 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -498,8 +498,8 @@ end
     A = CubicSpline(u, t; extrapolate = true)
 
     # Solution
-    P₁ = x -> 1 + 1.5x + x^2 + 0.5x^3 # for x ∈ [-1.0, 0.0]
-    P₂ = x -> 1 + 1.5x + x^2 - 0.5x^3 # for x ∈ [0.0, 1.0]
+    P₁ = x -> 1 + 1.5x + 0.75 * x^2 + 0.25 * x^3 # for x ∈ [-1.0, 0.0]
+    P₂ = x -> 1 + 1.5x + 0.75 * x^2 - 0.25 * x^3 # for x ∈ [0.0, 1.0]
 
     for (_t, _u) in zip(t, u)
         @test A(_t) == _u
@@ -534,8 +534,8 @@ end
     u = [0.0, 1.0, 3.0]
     t = [-1.0, 0.0, 1.0]
     A = CubicSpline(u, t; extrapolate = true)
-    @test A(-2.0) ≈ -2.0
-    @test A(2.0) ≈ 4.0
+    @test A(-2.0) ≈ -1.0
+    @test A(2.0) ≈ 5.0
     A = CubicSpline(u, t)
     @test_throws DataInterpolations.ExtrapolationError A(-2.0)
     @test_throws DataInterpolations.ExtrapolationError A(2.0)