diff --git a/examples/plot_narx.py b/examples/plot_narx.py
index d056962..f739528 100644
--- a/examples/plot_narx.py
+++ b/examples/plot_narx.py
@@ -44,22 +44,23 @@
 # %%
 # Build term libriary
 # -------------------
-#
 # To build a polynomial NARX model, it is normally have two steps:
 #
-# 1. Search the structure of the model, i.e., the terms in the model, e.g.,
-# :math:`u_0(k-1)u_0(k-3)`, :math:`u_0(k-2)u_1(k-3)`, etc.
-# 2. Learn the coefficients before the terms.
+# #. Search the structure of the model, i.e., the terms in the model, e.g.,
+#    :math:`u_0(k-1)u_0(k-3)`, :math:`u_0(k-2)u_1(k-3)`, etc.
+#
+# #. Learn the coefficients of the terms.
 #
 # To search the structure of the model, the candidate term libriary should be
-# constructed.
+# constructed by the following two steps.
+#
+# #. Time-shifted variables: the raw input-output data, i.e., :math:`u_0(k)`,
+#    :math:`u_1(k)`, and :math:`y(k)`, are converted into :math:`u_0(k-1)`,
+#    :math:`u_1(k-2)`, etc.
 #
-# 1. Time-shifted variables: the raw input-output data, i.e., :math:`u0(k)`,
-# :math:`u1(k)`, and :math:`y(k)`, are converted into :math:`u0(k-1)`, :math:`u1(k-2)`,
-# etc.
-# 2. Nonlinear terms: the time-shifted variables are onverted to nonlinear terms
-# via polynomial basis functions, e.g., :math:`u_0(k-1)^2`, :math:`u_0(k-1)u_0(k-3)`,
-# etc.
+# #. Nonlinear terms: the time-shifted variables are onverted to nonlinear terms
+#    via polynomial basis functions, e.g., :math:`u_0(k-1)^2`,
+#    :math:`u_0(k-1)u_0(k-3)`, etc.
 #
 #   .. rubric:: References
 #
@@ -116,6 +117,8 @@
 # ----------------
 # As the polynomical NARX is a linear function of the nonlinear tems,
 # the coefficient of each term can be easily estimated by oridnary least squares.
+# In the printed NARX model, it is found that :class:`FastCan` selects the correct
+# terms and the coefficients are close to the true values.
 
 from fastcan.narx import NARX, print_narx
 
@@ -126,9 +129,25 @@
 
 narx_model.fit(X, y)
 
-# In the printed NARX model, it is found that :class:`FastCan` selects the correct
-# terms and the coefficients are close to the true values.
 print_narx(narx_model)
+# %%
+# Automaticated NARX modelling workflow
+# -------------------------------------
+# We provide :meth:`narx.make_narx` to automaticate the workflow above.
+
+from fastcan.narx import make_narx
+
+auto_narx_model = make_narx(
+    X=X,
+    y=y,
+    n_features_to_select=4,
+    max_delay=3,
+    poly_degree=2,
+    verbose=0,
+).fit(X, y)
+
+print_narx(auto_narx_model)
+
 
 # %%
 # Plot NARX prediction performance