diff --git a/gbasis/boys.py b/gbasis/boys.py
new file mode 100644
index 00000000..09bc12de
--- /dev/null
+++ b/gbasis/boys.py
@@ -0,0 +1,64 @@
+"""Boys function for 2-electron integral evaluation.
+
+The Boys function F_m(x) = integral_0^1 t^(2m) * exp(-x * t^2) dt
+is the key auxiliary function in the Obara-Saika and Head-Gordon-Pople
+recurrence relations for 2-electron integrals.
+"""
+
+import numpy as np
+from scipy.special import hyp1f1, factorial2
+
+
+def boys_function(m, x):
+    """Evaluate the Boys function F_m(x).
+
+    The Boys function is defined as:
+        F_m(x) = integral_0^1 t^(2m) * exp(-x * t^2) dt
+
+    For small x, it is evaluated via the confluent hypergeometric function.
+    For large x, the asymptotic formula is used.
+
+    Parameters
+    ----------
+    m : int
+        Order of the Boys function (non-negative integer).
+    x : float or np.ndarray
+        Argument of the Boys function (non-negative).
+
+    Returns
+    -------
+    float or np.ndarray
+        Value of F_m(x).
+
+    Notes
+    -----
+    This function is central to the Obara-Saika (OS) recurrence relation
+    for evaluating 2-electron repulsion integrals [1]_.
+
+    References
+    ----------
+    .. [1] Obara, S.; Saika, A. J. Chem. Phys. 1986, 84, 3963.
+    """
+    x = np.asarray(x, dtype=float)
+    scalar = x.ndim == 0
+    x = np.atleast_1d(x)
+
+    result = np.empty_like(x)
+    large = x > 25.0  # threshold for asymptotic expansion
+
+    # Asymptotic formula for large x: F_m(x) ≈ (2m-1)!! / 2^(m+1) * sqrt(pi/x^(2m+1))
+    if np.any(large):
+        result[large] = (
+            np.math.factorial2(2 * m - 1)
+            / 2 ** (m + 1)
+            * np.sqrt(np.pi / x[large] ** (2 * m + 1))
+        )
+
+    # For small x: use F_m(x) = e^(-x) * sum or hypergeometric form
+    # F_m(x) = (1/(2m+1)) * 1F1(m+1/2; m+3/2; -x)
+    if np.any(~large):
+        result[~large] = (
+            1.0 / (2 * m + 1) * hyp1f1(m + 0.5, m + 1.5, -x[~large])
+        )
+
+    return result[0] if scalar else result
\ No newline at end of file
diff --git a/tests/test_boys.py b/tests/test_boys.py
new file mode 100644
index 00000000..ce99fc20
--- /dev/null
+++ b/tests/test_boys.py
@@ -0,0 +1,49 @@
+"""Tests for the Boys function utility."""
+
+import numpy as np
+import pytest
+from gbasis.boys import boys_function
+
+
+def test_boys_m0_x0():
+    """F_0(0) should equal 1."""
+    assert np.isclose(boys_function(0, 0.0), 1.0)
+
+
+def test_boys_m1_x0():
+    """F_1(0) should equal 1/3."""
+    assert np.isclose(boys_function(1, 0.0), 1.0 / 3.0)
+
+
+def test_boys_m2_x0():
+    """F_2(0) should equal 1/5."""
+    assert np.isclose(boys_function(2, 0.0), 1.0 / 5.0)
+
+
+def test_boys_m0_large_x():
+    """F_0(x) -> sqrt(pi) / (2*sqrt(x)) as x -> infinity."""
+    x = 100.0
+    expected = np.sqrt(np.pi) / (2.0 * np.sqrt(x))
+    assert np.isclose(boys_function(0, x), expected, rtol=1e-4)
+
+
+def test_boys_array_input():
+    """Boys function should handle array input and return correct shape."""
+    x = np.array([0.0, 1.0, 5.0, 30.0])
+    result = boys_function(0, x)
+    assert result.shape == (4,)
+
+
+def test_boys_positive_values():
+    """Boys function should return positive values for non-negative x."""
+    x = np.linspace(0.1, 50, 50)
+    for m in range(5):
+        assert np.all(boys_function(m, x) > 0)
+
+
+def test_boys_decreasing():
+    """F_m(x) should be strictly decreasing in x for fixed m."""
+    x = np.linspace(0.1, 10.0, 50)
+    for m in range(4):
+        vals = boys_function(m, x)
+        assert np.all(np.diff(vals) < 0)
\ No newline at end of file