diff --git a/.github/workflows/cookbook-practice.yml b/.github/workflows/cookbook-practice.yml
index a173140..906c266 100644
--- a/.github/workflows/cookbook-practice.yml
+++ b/.github/workflows/cookbook-practice.yml
@@ -27,7 +27,8 @@ jobs:
       run: |
         pip install --upgrade pip setuptools wheel
         pip install git+https://github.com/${{github.repository}}.git@${{ github.sha }}
+        pip install pytest-xdist
 
     - name: Run Tests
       run: |
-        pytest
+        pytest -n auto
diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml
index 0ac5c16..97c9c53 100644
--- a/.github/workflows/run-tests.yml
+++ b/.github/workflows/run-tests.yml
@@ -19,9 +19,7 @@ jobs:
         python-version: '3.11'
         cache: 'pip' # caching pip dependencies
     - run: |
-        pip install -r requirements.txt
-        cd docs/
-        pip install -r requirements.txt
-        cd ..
+        pip install -e .
+        pip install git+https://github.com/smallbunnies/myfunc.git@main
     - run: | # run both independent pytest and doctest
         pytest
diff --git a/Solverz/__init__.py b/Solverz/__init__.py
index 83e5e4d..fe26925 100644
--- a/Solverz/__init__.py
+++ b/Solverz/__init__.py
@@ -3,7 +3,7 @@
 from Solverz.equation.param import Param, IdxParam, TimeSeriesParam
 from Solverz.sym_algebra.symbols import idx, Para, iVar, iAliasVar
 from Solverz.sym_algebra.functions import (Sign, Abs, transpose, exp, Diag, Mat_Mul, sin, cos, Min, AntiWindUp,
-                                           Saturation, heaviside, ln)
+                                           Saturation, heaviside, ln, MulVarFunc, UniVarFunc)
 from Solverz.variable.variables import Vars, TimeVars, as_Vars
 from Solverz.solvers import *
 from Solverz.code_printer import made_numerical, module_printer
diff --git a/Solverz/code_printer/python/module/module_generator.py b/Solverz/code_printer/python/module/module_generator.py
index 120eebd..053d49b 100644
--- a/Solverz/code_printer/python/module/module_generator.py
+++ b/Solverz/code_printer/python/module/module_generator.py
@@ -67,8 +67,6 @@ def render_modules(eqs: SymEquations,
         code_dict["inner_Hvp"] = inner_Hvp['code_inner_Hvp']
         code_dict["sub_inner_Hvp"] = inner_Hvp['code_sub_inner_Hvp']
 
-
-
     def print_trigger_func_code():
         code_tfuc = dict()
         trigger_func = parse_trigger_func(eqs.PARAM)
@@ -176,7 +174,11 @@ def create_python_module(module_name,
 
 
 def print_init_code(eqn_type: str, module_name, eqn_param):
-    code = 'from .num_func import F_, J_\n'
+    code = '"""\n'
+    from ...._version import __version__
+    code += f'Python module generated by Solverz {__version__}\n'
+    code += '"""\n'
+    code += 'from .num_func import F_, J_\n'
     code += 'from .dependency import setting, p_, y\n'
     code += 'import time\n'
     match eqn_type:
@@ -207,7 +209,7 @@ def print_init_code(eqn_type: str, module_name, eqn_param):
     code += f'start = time.perf_counter()\n'
     code += code_compile.format(alpha=code_compile_args_F_J, beta=code_compile_args_Hvp)
     code += f'end = time.perf_counter()\n'
-    code += "print(f'Compiling time elapsed: {end-start}s')\n"
+    code += "print(f'Compiling time elapsed: {end - start}s')\n"
     return code
 
 
@@ -264,25 +266,14 @@ def print_module_code(code_dict: Dict[str, str], numba=False):
 
     return code
 
-#
-code_from_SolMuseum="""
-try:
-    import SolMuseum.num_api as SolMF
-except ImportError as e:
-    pass
 
-"""
 def print_dependency_code(modules):
     code = "import os\n"
     code += "current_module_dir = os.path.dirname(os.path.abspath(__file__))\n"
     code += 'from Solverz import load\n'
     code += 'auxiliary = load(f"{current_module_dir}\\\\param_and_setting.pkl")\n'
     code += 'from numpy import *\n'
-    code += 'import numpy as np\n'
-    code += 'import Solverz.num_api.custom_function as SolCF\n'  # import Solverz built-in func
-    code += code_from_SolMuseum
-    code += 'import scipy.sparse as sps\n'
-    code += 'from numba import njit\n'
+    code += 'from Solverz.num_api.module_parser import *\n'
     code += 'setting = auxiliary["eqn_param"]\n'
     code += 'row = setting["row"]\n'
     code += 'col = setting["col"]\n'
diff --git a/Solverz/code_printer/python/module/test/test_module_generator.py b/Solverz/code_printer/python/module/test/test_module_generator.py
index 0731404..1de3030 100644
--- a/Solverz/code_printer/python/module/test/test_module_generator.py
+++ b/Solverz/code_printer/python/module/test/test_module_generator.py
@@ -18,16 +18,7 @@
 from Solverz import load
 auxiliary = load(f"{current_module_dir}\\param_and_setting.pkl")
 from numpy import *
-import numpy as np
-import Solverz.num_api.custom_function as SolCF
-
-try:
-    import SolMuseum.num_api as SolMF
-except ImportError as e:
-    pass
-
-import scipy.sparse as sps
-from numba import njit
+from Solverz.num_api.module_parser import *
 setting = auxiliary["eqn_param"]
 row = setting["row"]
 col = setting["col"]
@@ -38,7 +29,12 @@
 y = auxiliary["vars"]
 """
 
-expected_init = r"""from .num_func import F_, J_
+from ....._version import __version__
+
+expected_init = r'''"""
+Python module generated by Solverz {vs}
+"""
+from .num_func import F_, J_
 from .dependency import setting, p_, y
 import time
 from Solverz.num_api.num_eqn import nAE
@@ -60,8 +56,8 @@
     v = ones_like(y)
     mdl.HVP(y, p_, v)
 end = time.perf_counter()
-print(f'Compiling time elapsed: {end-start}s')
-"""
+print(f'Compiling time elapsed: {{end - start}}s')
+'''.format(vs=__version__)
 
 
 def test_AE_module_printer():
diff --git a/Solverz/num_api/module_parser.py b/Solverz/num_api/module_parser.py
index 7c0fd44..513f3d0 100644
--- a/Solverz/num_api/module_parser.py
+++ b/Solverz/num_api/module_parser.py
@@ -4,13 +4,28 @@
 import warnings
 
 import Solverz.num_api.custom_function as SolCF
+import numpy as np
+import scipy.sparse as sps
+from numba import njit
 
-modules = [{'SolCF': SolCF, 'np': numpy, 'sps': scipy.sparse}, 'numpy']
-# We preserve the 'numpy' here in case one uses functions from sympy instead of from Solverz
+module_dict = {'SolCF': SolCF, 'np': np, 'sps': sps, 'njit': njit}
 
 # parse modules from museum
 try:
     import SolMuseum.num_api as SolMF
-    modules[0]['SolMF'] = SolMF
+    module_dict['SolMF'] = SolMF
 except ModuleNotFoundError as e:
     warnings.warn(f'Failed to import num api from SolMuseum: {e}')
+
+# parse user defined functions
+
+try:
+    import myfunc
+    print('User module detected.')
+    module_dict['myfunc'] = myfunc
+except ModuleNotFoundError as e:
+    pass
+
+modules = [module_dict, 'numpy']
+# We preserve the 'numpy' here in case one uses functions from sympy instead of from Solverz
+__all__ = list(module_dict.keys())
diff --git a/Solverz/num_api/test/test_udm.py b/Solverz/num_api/test/test_udm.py
new file mode 100644
index 0000000..3cb00a7
--- /dev/null
+++ b/Solverz/num_api/test/test_udm.py
@@ -0,0 +1,47 @@
+"""
+Test the user defined modules.
+"""
+
+
+def test_udm():
+
+    from Solverz import Model, Var, Eqn, made_numerical, MulVarFunc
+    import numpy as np
+
+    class Min(MulVarFunc):
+        arglength = 2
+
+        def fdiff(self, argindex=1):
+            if argindex == 1:
+                return dMindx(*self.args)
+            elif argindex == 2:
+                return dMindy(*self.args)
+
+        def _numpycode(self, printer, **kwargs):
+            return (f'myfunc.Min' + r'(' +
+                    ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
+
+    class dMindx(MulVarFunc):
+        arglength = 2
+
+        def _numpycode(self, printer, **kwargs):
+            return (f'myfunc.dMindx' + r'(' +
+                    ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
+
+    class dMindy(MulVarFunc):
+        arglength = 2
+
+        def _numpycode(self, printer, **kwargs):
+            return (f'myfunc.dMindy' + r'(' +
+                    ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
+
+    m = Model()
+    m.x = Var('x', [1, 2])
+    m.y = Var('y', [3, 4])
+    m.f = Eqn('f', Min(m.x, m.y))
+    sae, y0 = m.create_instance()
+    ae = made_numerical(sae, y0, sparse=True)
+    np.testing.assert_allclose(ae.F(y0, ae.p), np.array([1.0, 2.0]))
+    np.testing.assert_allclose(ae.J(y0, ae.p).toarray(), np.array([[1., 0., 0., 0.],
+                                                                   [0., 1., 0., 0.]]))
+
diff --git a/docs/src/advanced.md b/docs/src/advanced.md
index 95e8227..5d795d6 100644
--- a/docs/src/advanced.md
+++ b/docs/src/advanced.md
@@ -4,131 +4,198 @@
 
 ## Writing Custom Functions
 Sometimes one may have functions that go beyond the Solverz built-in library. This guide will describe how to create 
-such custom functions in Solverz, so that the functions can be incorporated into numerical simulations. The philosophy 
-of function customization comes from Sympy, it helps to learn the [Sympy basics](https://docs.sympy.org/latest/index.html) 
+such custom functions using Solverz and inform Solverz of their paths, so that the functions can be incorporated into 
+numerical simulations. 
+
+```{note}
+Alternatively, one can directly contribute to the [SolMuseum](https://solmuseum.solverz.org/stable/) 
+library so that 1) others can utilize your models/functions and 2) one avoids configuring the module paths.
+```
+
+```{note}
+The philosophy of function customization comes from Sympy, it helps to learn the [Sympy basics](https://docs.sympy.org/latest/index.html) 
 and read the [Sympy tutorial of custom functions](https://docs.sympy.org/latest/guides/custom-functions.html) for an overview.
+```
 
-As a motivating example for this document, let's create a custom function class representing the $\min$ function. We 
-want use $\min$ to determine the smaller one of two operands, which can be defined by
+```{note}
+In Solverz, the numerical computations are mainly dependent on the prevailing numerical libraries such as numpy and scipy. 
+It is recommended that one first gets familiar with the [numpy](https://numpy.org/doc/stable/user/index.html) and 
+[scipy](https://docs.scipy.org/doc/scipy/tutorial/index.html).
+```
+
+### An Illustrative Example
+
+As a motivating example for this document, let's create a custom function class representing the $\min$ function.
+The $\min$ function is typical in controllers of many industrial applications, which can be defined by
 
 ```{math}
 \min(x,y)=
 \begin{cases}
 x&x\leq y\\
 y& x>y
-\end{cases}
+\end{cases}.
 ```
 
-We also want to extend the function to vector input, that is, capable of finding the element-wise minimum of two vectors.
-To summarize, we shall implement $\min$ that
+To incorporate $\min$ in our simulation modelling, its symbolic and numerical implementations shall be defined. 
+Specifically,
 
-1. evaluates $\min(x,y)$ correctly
-2. can be derived proper derivatives with respect to $x$ and $y$.
+1. a symbolic function `min` can be called to represent the $\min$ function;
+2. the symbolic derivatives of `min` are automatically derived for the Jacobian block parser;
+3. the numerical interfaces are defined so that the `min` function and its derivatives can be correctly evaluated.
 
-However, it is difficult to devise the analytical derivatives of $\min$. We should perform the trick that rewrites 
-$\min(x,y)$ as
+First, we define the numerical interfaces. The derivatives of $\min$ function are
 
 ```{math}
-x*\operatorname{lessthan}(x,y)+y*(1-\operatorname{lessthan}(x,y)).
+\pdv{\min(x,y)}{x}=
+\begin{cases}
+1&x\leq y\\
+0& x>y
+\end{cases},\quad
+\pdv{\min(x,y)}{y}=
+\begin{cases}
+0&x\leq y\\
+1& x>y
+\end{cases}.
 ```
 
-where the $\operatorname{lessthan}(x,y)$ function mathematically denotes the $\leq$ operator and returns 1 if 
-$x\leq y$ else 0. Since $\operatorname{lessthan}(x,y)$ can only be either 1 or 0, the above transformation holds. 
+Let us create a `myfunc` directory and put the numerical codes in the `myfunc.py` file that looks like
 
-If the derivatives of $\operatorname{lessthan}(x,y)$ with respect to any argument are zero, then we have
-```{math}
-\frac{\partial}{\partial x}\min(x,y)=
-\operatorname{lessthan}(x,y)
+```python
+# myfunc.py
+import numpy as np
+from numba import njit
+
+
+@njit(cache=True)
+def Min(x, y):
+    x = np.asarray(x).reshape((-1,))
+    y = np.asarray(y).reshape((-1,))
+
+    z = np.zeros_like(x)
+
+    for i in range(len(x)):
+        if x[i] <= y[i]:
+            z[i] = x[i]
+        else:
+            z[i] = y[i]
+    return z
+
+
+@njit(cache=True)
+def dMindx(x, y):
+    x = np.asarray(x).reshape((-1,))
+    y = np.asarray(y).reshape((-1,))
+
+    z = np.zeros_like(x)
+
+    for i in range(len(x)):
+        if x[i] <= y[i]:
+            z[i] = 1
+        else:
+            z[i] = 0
+    return z
+
+
+@njit(cache=True)
+def dMindy(x, y):
+    return 1-dMindx(x, y)
 ```
-and
-```{math}
-\frac{\partial}{\partial y}\min(x,y)=
-1-\operatorname{lessthan}(x,y).
+
+In `myfunc.py`, we use `Min` to avoid conflicting with the built-in `min` function. 
+The `@njit(cache)` decorator is used to perform the jit-compilation and hence speed up the numerical codes.
+
+Then let us install the `myfunc` module, so that Solverz can import the `myfunc` module. Use the terminal to switch
+to the `myfunc` module directory. Add a `pyproject.toml` file there. 
+
+```{note}
+One can clone the `pyproject.toml` file from [example repo](https://github.com/smallbunnies/myfunc).
+```
+
+Use the following command to install the module.
+
+```shell
+pip install -e .
+```
+
+```{note}
+Because `myfunc` module is installed in the editable mode, one can change the numerical implementations in `myfunc.py`
+with great freedom.
 ```
-Hence, it suffices to have a custom $\operatorname{lessthan}(x,y)$ function that
 
-1. evaluates $\operatorname{lessthan}(x,y)$ correctly
-2. has zero-derivative with respect to $x$ or $y$.
+As for the symbolic implementation, let us start by creating a `Min.py` file and subclassing `MulVarFunc` there with
 
-Let us start by subclassing `MulVarFunc`
 ```python
-from Solverz.sym_algebra.functions import MulVarFunc
+from Solverz import MulVarFunc
+
 class Min(MulVarFunc):
     pass
-class LessThan(MulVarFunc):
+
+class dMindx(MulVarFunc):
+    pass
+
+class dMindy(MulVarFunc):
     pass
 ```
-The `MulVarFunc` is the base class of multi-variate functions in Solverz. 
-At this point, `Min` has no behaviors defined on it. To automatically evaluate the `Min` function, we ought to define 
-the **_class method_** `eval()`. `eval()` should take the arguments of the function and return the value 
-$x*\operatorname{lessthan}(x,y)+y*(1-\operatorname{lessthan}(x,y))$:
+
+The `MulVarFunc` is the base class of symbolic multi-variate functions in Solverz. 
+
+At this point, `Min` and its derivatives have no behaviors defined on it. To instruct Solverz in the differentiation 
+rule of `Min` and the numerical implementations, we shall add following lines
 ```python
 class Min(MulVarFunc):
-    @classmethod
-    def eval(cls, x, y):
-        return x * LessThan(x, y) + y * (1 - LessThan(x, y))
-```
-```python
->>> from Solverz import Var
->>> Min(Var('x',0),Var('y',0))
-... x*((x)<=(y)) + y*(1 - ((x)<=(y)))
-```
-To define the differentiation of  `LessThan()`, we have
-```python
-from sympy import Integer
-class LessThan(MulVarFunc):
-    """
-    Represents < operator
-    """
+    arglength = 2
 
-    def _eval_derivative(self, s):
-        return Integer(0)
+    def fdiff(self, argindex=1):
+        if argindex == 1:
+            return dMindx(*self.args)
+        elif argindex == 2:
+            return dMindy(*self.args)
+
+    def _numpycode(self, printer, **kwargs):
+        return (f'myfunc.Min' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-    def _sympystr(self, printer, **kwargs):
-        return '(({op1})<=({op2}))'.format(op1=printer._print(self.args[0]),
-                                           op2=printer._print(self.args[1]))
+
+class dMindx(MulVarFunc):
+    arglength = 2
 
     def _numpycode(self, printer, **kwargs):
-        return r'SolLessThan(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
+        return (f'myfunc.dMindx' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-    def _lambdacode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
 
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-```
-Here, the `_sympystr()` method defines its string representation:
-```python
->>> LessThan(Var('x'), Var('y'))
-... ((x)<=(y))
-```
-The `_eval_derivative()` method forces the derivatives of `LessThan()` to be zero:
-```python
-from Solverz import iVar
->>> Min(Var('x',0),Var('y',0)).diff(iVar('x'))
-... ((x)<=(y))
-```
-where `iVar` is the internal variable type of Solverz, `diff()` is the method to derive derivatives.
+class dMindy(MulVarFunc):
+    arglength = 2
+
+    def _numpycode(self, printer, **kwargs):
+        return (f'myfunc.dMindy' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
-The `_numpycode()` function defines what should `LessThan()` be printed to  in numerical codes. Here, we define the 
-`SolLessThan()` as the numerical implementation of  `LessThan()`. Given array `[0,2,-1]` and `[1,2,3]`:
-```python
->>> import numpy as np
->>> SolLessThan(np.array([0, 2, -1]), np.array([1,2,3]))
-... array([1, 0, 1])
-```
-```{note}
-In Solverz, the numerical computations are mainly dependent on the prevailing numerical libraries such as numpy and scipy. 
-It is recommended that one first gets familiar with the [numpy](https://numpy.org/doc/stable/user/index.html) and 
-[scipy](https://docs.scipy.org/doc/scipy/tutorial/index.html).
 ```
-The implementation of `SolLessThan()` should be put in the `Solverz.num_api.custom_function` module:
+
+where the `fdiff` function should return the derivative of the function, without considering the chain rule, 
+with respect to the `argindex`-th variable; the `_numpycode` functions define the numerical implementations of the 
+functions. Since the `myfunc` module has been installed, the numerical implementations can be called by
+`myfunc.Min`.
+
+After finish the above procedures, we can finally use the `Min` function in our simulation modelling. An example is as 
+follows.
+
 ```python
-@implements_nfunc('SolLessThan')
-@njit(cache=True)
-def SolLessThan(x, y):
-    x = np.asarray(x).reshape((-1,))
-    return (x < y).astype(np.int32)
+from Solverz import Model, Var, Eqn, made_numerical
+from Min import Min
+
+m = Model()
+m.x = Var('x', [1, 2])
+m.y = Var('y', [3, 4])
+m.f = Eqn('f', Min(m.x, m.y))
+sae, y0 = m.create_instance()
+ae = made_numerical(sae, y0, sparse=True)
 ```
-The `implements_nfunc()` cannot be omitted and the `njit()` decorator enables the numba-based dynamic compilation for efficiency.
 
+We will have the output
+
+```shell
+>>> ae.F(y0, ae.p)
+array([1.0, 2.0])
+```