diff --git a/docs/code/sensitivity/chatterjee/local_exponential.py b/docs/code/sensitivity/chatterjee/local_exponential.py
new file mode 100644
index 00000000..1fd0ef0d
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/local_exponential.py
@@ -0,0 +1,20 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+
+
+def evaluate(X: np.array) -> np.array:
+    r"""A non-linear function that is used to demonstrate sensitivity index.
+
+    .. math::
+        f(x) = \exp(x_1 + 2*x_2)
+    """
+
+    Y = np.exp(X[:, 0] + 2 * X[:, 1])
+
+    return Y
diff --git a/docs/code/sensitivity/chatterjee/local_ishigami.py b/docs/code/sensitivity/chatterjee/local_ishigami.py
new file mode 100644
index 00000000..e5af649f
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/local_ishigami.py
@@ -0,0 +1,23 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+
+
+def evaluate(X, params=[7, 0.1]):
+    """Non-monotonic Ishigami-Homma three parameter test function"""
+
+    a = params[0]
+    b = params[1]
+
+    Y = (
+        np.sin(X[:, 0])
+        + a * np.power(np.sin(X[:, 1]), 2)
+        + b * np.power(X[:, 2], 4) * np.sin(X[:, 0])
+    )
+
+    return Y
diff --git a/docs/code/sensitivity/chatterjee/local_sobol_func.py b/docs/code/sensitivity/chatterjee/local_sobol_func.py
new file mode 100644
index 00000000..1ccabc6d
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/local_sobol_func.py
@@ -0,0 +1,42 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+
+import numpy as np
+import copy
+
+
+def evaluate(X, a_values):
+
+    dims = len(a_values)
+    g = 1
+
+    for i in range(dims):
+        g_i = (np.abs(4 * X[:, i] - 2) + a_values[i]) / (1 + a_values[i])
+        g *= g_i
+
+    return g
+
+
+def sensitivities(a_values):
+
+    dims = len(a_values)
+
+    Total_order = np.zeros((dims, 1))
+
+    V_i = (3 * (1 + a_values) ** 2) ** (-1)
+
+    total_variance = np.prod(1 + V_i) - 1
+
+    First_order = V_i / total_variance
+
+    for i in range(dims):
+
+        rem_First_order = copy.deepcopy(V_i)
+        rem_First_order[i] = 0
+        Total_order[i] = V_i[i] * np.prod(rem_First_order + 1) / total_variance
+
+    return First_order.reshape(-1, 1), Total_order
diff --git a/docs/code/sensitivity/chatterjee/plot_chatterjee_exponential.py b/docs/code/sensitivity/chatterjee/plot_chatterjee_exponential.py
new file mode 100644
index 00000000..2922b97a
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/plot_chatterjee_exponential.py
@@ -0,0 +1,54 @@
+"""
+
+Exponential function
+==============================================
+
+.. math::
+    f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
+
+"""
+
+# %%
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Normal
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.chatterjee import Chatterjee
+
+# %%
+# Create Model object
+model = PythonModel(
+    model_script="local_exponential.py",
+    model_object_name="evaluate",
+    var_names=[
+        "X_1",
+        "X_2",
+    ],
+    delete_files=True,
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Normal(0, 1)] * 2)
+
+# %% [markdown]
+# Compute Chatterjee indices
+
+# %%
+SA = Chatterjee(runmodel_obj, dist_object)
+
+# Compute Sobol indices using the pick and freeze algorithm
+computed_indices = SA.run(n_samples=1_000_000)
+
+# %% [markdown]
+# Cramer-von Mises sensitivity analysis
+#
+# Expected value of the sensitivity indices:
+#
+# $S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145$
+#
+# $S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693$
+
+# %%
+computed_indices["chatterjee_i"]
diff --git a/docs/code/sensitivity/chatterjee/plot_chatterjee_ishigami.py b/docs/code/sensitivity/chatterjee/plot_chatterjee_ishigami.py
new file mode 100644
index 00000000..66897d67
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/plot_chatterjee_ishigami.py
@@ -0,0 +1,58 @@
+r"""
+
+Ishigami function
+==============================================
+
+.. math::
+    f(x_1, x_2, x_3) = sin(x_1) + a \cdot sin^2(x_2) + b \cdot x_3^4 sin(x_1)
+
+.. math::
+    x_1, x_2, x_3 \sim \mathcal{U}(-\pi, \pi), \quad a, b\in \mathbb{R}
+
+"""
+
+# %%
+import numpy as np
+
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Uniform
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.chatterjee import Chatterjee
+
+# %%
+# Create Model object
+model = PythonModel(
+    model_script="local_ishigami.py",
+    model_object_name="evaluate",
+    var_names=[r"$X_1$", "$X_2$", "$X_3$"],
+    delete_files=True,
+    params=[7, 0.1],
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Uniform(-np.pi, 2 * np.pi)] * 3)
+
+# %% [markdown]
+# Compute Chatterjee indices
+
+# %%
+SA = Chatterjee(runmodel_obj, dist_object)
+
+computed_indices = SA.run(
+    n_samples=100_000,
+    estimate_sobol_indices=True,
+    num_bootstrap_samples=100,
+    confidence_level=0.95,
+)
+
+# %%
+computed_indices["chatterjee_i"]
+
+# %%
+computed_indices["CI_chatterjee_i"]
+
+# %%
+computed_indices["sobol_i"]
diff --git a/docs/code/sensitivity/chatterjee/plot_chatterjee_sobol_func.py b/docs/code/sensitivity/chatterjee/plot_chatterjee_sobol_func.py
new file mode 100644
index 00000000..57813142
--- /dev/null
+++ b/docs/code/sensitivity/chatterjee/plot_chatterjee_sobol_func.py
@@ -0,0 +1,70 @@
+r"""
+
+Sobol function
+==============================================
+
+.. math::
+
+    g(x_1, x_2, \ldots, x_D) := \prod_{i=1}^{D} \frac{|4x_i - 2| + a_i}{1 + a_i},
+
+where,
+
+.. math::
+    x_i \sim \mathcal{U}(0, 1), \quad a_i \in \mathbb{R}.
+
+"""
+
+# %%
+import numpy as np
+
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+from UQpy.distributions import Uniform
+from UQpy.distributions.collection.JointIndependent import JointIndependent
+from UQpy.sensitivity.chatterjee import Chatterjee
+
+# %%
+# Create Model object
+num_vars = 6
+a_vals = np.array([0.0, 0.5, 3.0, 9.0, 99.0, 99.0])
+
+model = PythonModel(
+    model_script="local_sobol_func.py",
+    model_object_name="evaluate",
+    var_names=["X_" + str(i) for i in range(num_vars)],
+    delete_files=True,
+    a_values=a_vals,
+)
+
+runmodel_obj = RunModel(model=model)
+
+# Define distribution object
+dist_object = JointIndependent([Uniform(0, 1)] * num_vars)
+
+# %% [markdown]
+# Compute Chatterjee indices
+
+# %%
+SA = Chatterjee(runmodel_obj, dist_object)
+
+# Compute Sobol indices using the pick and freeze algorithm
+computed_indices = SA.run(n_samples=500_000, estimate_sobol_indices=True)
+
+# %%
+computed_indices["chatterjee_i"]
+
+# %% [markdown]
+# $S_1$ = 5.86781190e-01
+#
+# $S_2$ = 2.60791640e-01
+#
+# $S_3$ = 3.66738244e-02
+#
+# $S_4$ = 5.86781190e-03
+#
+# $S_5$ = 5.86781190e-05
+#
+# $S_6$ = 5.86781190e-05
+
+# %%
+computed_indices["sobol_i"]