diff --git a/docs/src/interface.md b/docs/src/interface.md
index b7c51337..f9b558fd 100644
--- a/docs/src/interface.md
+++ b/docs/src/interface.md
@@ -32,11 +32,64 @@ A2(300.0)
 ```
 
 !!! note
-    
+
     The values computed beyond the range of the time points provided during interpolation will not be reliable, as these methods only perform well within the range and the first/last piece polynomial fit is extrapolated on either side which might not reflect the true nature of the data.
 
 The keyword `cache_parameters = true` can be passed to precalculate parameters at initialization, making evaluations cheaper to compute. This is not compatible with modifying `u` and `t`. The default `cache_parameters = false` does however not prevent allocation in every interpolation constructor call.
 
+### In-place Evaluation (Allocation-free)
+
+When performance is critical, such as in ODE solvers or tight loops, you can use the in-place variant to avoid memory allocations. This is particularly useful when you need to interpolate at many points repeatedly.
+
+To use in-place interpolation, pass a pre-allocated output array as the first argument:
+
+```@example interface
+# Pre-allocate output array
+t_eval = [50.0, 100.0, 150.0, 200.0]
+u_out = similar(u, length(t_eval))
+
+# In-place interpolation: interp(output, t_values)
+A1(u_out, t_eval)
+
+u_out
+```
+
+The in-place form `interp(out, t)` writes the interpolated values directly into `out`, avoiding allocation of a new array. The output array must have the same length as the input time vector.
+
+For vector-of-arrays data (where `u` is a `Vector{Vector}` or `Vector{Matrix}`), the output should be a pre-allocated vector of arrays with the same structure:
+
+```@example interface
+# Vector of vectors example
+u_vov = [[14.7, 7.35], [11.51, 5.76], [10.41, 5.21], [14.95, 7.48], [12.24, 6.12], [11.22, 5.61]]
+A_vov = CubicSpline(u_vov, t)
+
+# Pre-allocate output as Vector{Vector}
+t_eval_3 = [50.0, 100.0, 150.0]
+out_vov = [zeros(2) for _ in 1:length(t_eval_3)]
+A_vov(out_vov, t_eval_3)
+
+out_vov
+```
+
+For multi-dimensional data (where `u` is a matrix or higher-dimensional array), the output array should have the same leading dimensions as `u`, with the last dimension matching the length of `t`:
+
+```@example interface
+# Matrix example (stacked form)
+u_matrix = [14.7 11.51 10.41 14.95 12.24 11.22;
+            7.35 5.76 5.21 7.48 6.12 5.61]
+A_matrix = CubicSpline(u_matrix, t)
+
+# Pre-allocate for 3 evaluation points
+out_matrix = zeros(2, 3)
+A_matrix(out_matrix, t_eval_3)
+
+out_matrix
+```
+
+!!! tip "Performance tip"
+
+    Using in-place interpolation can significantly reduce allocations in performance-critical code. For example, in an ODE derivative function, switching from `interp(t_vec)` to `interp(out, t_vec)` can eliminate allocations entirely within the hot loop.
+
 ## Derivatives
 
 Derivatives of the interpolated curves can also be computed at any point for all the methods. Derivatives upto second order is supported where first order derivative is computed analytically and second order using `ForwardDiff.jl`. Order is passed as the third argument. It is 1 by default.