diff --git a/deriv/README.md b/deriv/README.md
index 096613e8d..b66312a63 100644
--- a/deriv/README.md
+++ b/deriv/README.md
@@ -63,6 +63,18 @@ for values greater than `x`. 📉
 
 This function is equivalent to calling `deriv.forward` with a negative step-size.
 
+### `partial`
+
+```v ignore
+pub fn partial(f fn ([]f64) f64, x []f64, variable int, h f64) (f64, f64)
+```
+
+This function computes the partial derivative of the function `f` with respect to
+a specified variable at point `x` using step-size `h`. It returns the derivative
+in `result` and an error estimate in `abserr`. The function `f` should take an array
+of coordinates and return a single value. This method provides both the derivative
+and its error estimate.
+
 ## References and Further Reading
 
 This work is a spiritual descendent of the Differentiation module in [GSL](https://github.com/ampl/gsl). 📖
diff --git a/deriv/deriv.v b/deriv/deriv.v
index 47a2b33cd..eb53c0758 100644
--- a/deriv/deriv.v
+++ b/deriv/deriv.v
@@ -1,6 +1,7 @@
 module deriv
 
 import vsl.func
+import vsl.errors
 import vsl.internal.prec
 import math
 
@@ -113,3 +114,32 @@ pub fn forward(f func.Fn, x f64, h f64) (f64, f64) {
 pub fn backward(f func.Fn, x f64, h f64) (f64, f64) {
 	return forward(f, x, -h)
 }
+
+/**
+* partial is a function that computes the partial derivative of a multivariable function with respect to a specified variable.
+*
+* @param f       The multivariable function for which the partial derivative is to be computed.
+* @param x       The point at which the partial derivative is to be computed, represented as an array of coordinates.
+* @param variable The index of the variable with respect to which the partial derivative is to be computed.
+* @param h       The step size to be used in the central difference method.
+*
+* @return        A tuple containing the value of the partial derivative and the estimated error.
+*/
+pub fn partial(f fn ([]f64) f64, x []f64, variable int, h f64) (f64, f64) {
+	if variable < 0 || variable >= x.len {
+		errors.vsl_panic('Invalid variable index', .efailed)
+	}
+
+	// Define a helper function that converts the multivariate function
+	// to a univariate function for the specified variable
+	partial_helper := func.Fn{
+		f: fn [f, x, variable] (t f64, _ []f64) f64 {
+			mut x_new := x.clone()
+			x_new[variable] = t
+			return f(x_new)
+		}
+	}
+
+	// Use the central difference method to compute the partial derivative
+	return central(partial_helper, x[variable], h)
+}
diff --git a/deriv/deriv_test.v b/deriv/deriv_test.v
index f25a5573f..041c53d79 100644
--- a/deriv/deriv_test.v
+++ b/deriv/deriv_test.v
@@ -68,6 +68,18 @@ fn df6(x f64, _ []f64) f64 {
 	return -1.0 / (x * x)
 }
 
+fn f_multi(x []f64) f64 {
+	return x[0] * x[0] + x[1] * x[1] // f(x,y) = x^2 + y^2
+}
+
+fn df_multi_dx(x []f64) f64 {
+	return 2 * x[0] // ∂f/∂x = 2x
+}
+
+fn df_multi_dy(x []f64) f64 {
+	return 2 * x[1] // ∂f/∂y = 2y
+}
+
 fn test_deriv() {
 	f1_ := func.Fn.new(f: f1)
 	df1_ := func.Fn.new(f: df1)
@@ -100,6 +112,18 @@ fn test_deriv() {
 	assert deriv_test('central', f6_, df6_, 10.0)
 	assert deriv_test('forward', f6_, df6_, 10.0)
 	assert deriv_test('backward', f6_, df6_, 10.0)
+
+	// Partial derivative test
+	x := [2.0, 3.0]
+	h := 1e-5
+
+	// Partial derivative with respect to x
+	dx, _ := partial(f_multi, x, 0, h)
+	assert float64.tolerance(dx, df_multi_dx(x), 1e-5)
+
+	// Partial derivative with respect to y
+	dy, _ := partial(f_multi, x, 1, h)
+	assert float64.tolerance(dy, df_multi_dy(x), 1e-5)
 }
 
 fn deriv_test(deriv_method string, f func.Fn, df func.Fn, x f64) bool {