diff --git a/docs/examples/emulator/robustgasp/plot_ppgasp.py b/docs/examples/emulator/robustgasp/plot_ppgasp.py
new file mode 100644
index 0000000000000000000000000000000000000000..517dd880a6e4ef68aa6d92431c35aafe6f16eb04
--- /dev/null
+++ b/docs/examples/emulator/robustgasp/plot_ppgasp.py
@@ -0,0 +1,109 @@
+"""
+
+PPGaSP: GP emulation for multi-output functions
+===============================================
+"""
+
+
+# %% md
+# 
+# This example shows how to apply Gaussian process emulation to a multi-output
+# function using class :class:`.PPGaSP`.
+#
+# The multi-output function that we are going to look at is the `DIAMOND`
+# (diplomatic and military operations in a non-warfighting domain) computer model.
+# It is used as a testbed to illustrate the `PP GaSP` emulator in the `R` package
+# `RobustGaSP`, see :cite:t:`Gu2019` for more detail.
+#
+# The simulator has :math:`13` input parameters and :math:`5` outputs. Namely,
+# :math:`\mathbf{y}=f(\mathbf{x})` where :math:`\mathbf{x}=(x_1,\ldots,x_{13})^T`
+# and :math:`\mathbf{y}=(y_1,\ldots,y_5)^T`.
+#
+# The training and testing data are provided in the folder '.../tests/data/'.
+# We first load the training and testing data. You may need to modify ``dir_data```
+# according to where you save them.
+
+import numpy as np
+import os
+
+dir_data = os.path.abspath('../../../../tests/data/')
+
+humanityX = np.genfromtxt(os.path.join(dir_data, 'humanityX.csv'), delimiter=',')
+humanityY = np.genfromtxt(os.path.join(dir_data, 'humanityY.csv'), delimiter=',')
+print(f"Number of training data points: ", humanityX.shape[0])
+print(f"Input dimension: ", humanityX.shape[1])
+print(f"Output dimension: ", humanityY.shape[1])
+
+humanityXt = np.genfromtxt(os.path.join(dir_data, 'humanityXt.csv'), delimiter=',')
+humanityYt = np.genfromtxt(os.path.join(dir_data, 'humanityYt.csv'), delimiter=',')
+print(f"Number of testing data points: ", humanityXt.shape[0])
+
+# %% md
+# ``humanityX`` and ``humanitY`` are the training data, corresponding to ``design``
+# and ``response`` respectively. ``humanityXt`` are testing input data, at which
+# we are going to make predictions once the emulator is trained. ``humanityYt``
+# are the true outputs at ``humanityXt``, which is then used to validate the
+# performance of the trained emulator.
+
+# %% md
+# 
+# To build a `PP GaSP` emulator for the above simulator,  first import class
+# :class:`.PPGaSP` by
+
+from psimpy.emulator import PPGaSP
+
+# %% md
+# 
+# Then, create an instance of :class:`.PPGaSP`.  The parameter ``ndim``
+# (dimension of function input ``x``) must be specified. Optional parameters, such
+# as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave
+# all the optional parameters to their default values.
+
+emulator = PPGaSP(ndim=humanityX.shape[1])
+
+
+# %% md
+# 
+# Next, we train the `PP GaSP` emulator based on the training data using
+# :py:meth:`.PPGaSP.train` method.
+
+emulator.train(design=humanityX, response=humanityY)
+
+# %% md
+#
+# With the trained emulator, we can make predictions for any arbitrary set of
+# input points using :py:meth:`PPGaSP.predict` method.
+# Here, we make predictions at testing input points ``humanityXt``.
+
+predictions = emulator.predict(humanityXt)
+
+# %% md
+# 
+# We can validate the performance of the trained emulator based on the true outputs
+# ``humanityYt``.
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(5, 1, figsize=(6,15))
+
+for i in range(humanityY.shape[1]):
+
+    ax[i].set_xlabel(f'Actual $y_{i+1}$')
+    ax[i].set_ylabel(f'Emulator-predicted $y_{i+1}$')
+    ax[i].set_xlim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))
+    ax[i].set_ylim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))
+
+    _ = ax[i].plot([np.min(humanityYt[:,i]),np.max(humanityYt[:,i])], [np.min(humanityYt[:,i]),np.max(humanityYt[:,i])])
+    _ = ax[i].errorbar(humanityYt[:,i], predictions[:,i,0], predictions[:,i,3], fmt='.', linestyle='', label='prediction and std')
+    _ = ax[i].legend()
+
+plt.tight_layout()
+
+
+# %%md
+#
+# We can also draw any number of samples at testing input ``humanityXt`` using 
+# :py:meth:`.PPGaSPsample()` method.
+
+samples = emulator.sample(humanityXt, nsamples=10)
+print("Shape of samples: ", samples.shape)
\ No newline at end of file