diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml
index 60488e1..0ac5c16 100644
--- a/.github/workflows/run-tests.yml
+++ b/.github/workflows/run-tests.yml
@@ -24,4 +24,4 @@ jobs:
         pip install -r requirements.txt
         cd ..
     - run: | # run both independent pytest and doctest
-        pytest -vv
+        pytest
diff --git a/Solverz/__init__.py b/Solverz/__init__.py
index 70b1125..83e5e4d 100644
--- a/Solverz/__init__.py
+++ b/Solverz/__init__.py
@@ -1,10 +1,9 @@
-from Solverz.equation.eqn import Eqn, Ode, HyperbolicPde
+from Solverz.equation.eqn import Eqn, Ode
 from Solverz.equation.equations import AE, FDAE, DAE
 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)
-from Solverz.num_api.custom_function import minmod_flag, minmod
 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/inline/inline_printer.py b/Solverz/code_printer/python/inline/inline_printer.py
index 0fd2663..b66fddc 100644
--- a/Solverz/code_printer/python/inline/inline_printer.py
+++ b/Solverz/code_printer/python/inline/inline_printer.py
@@ -4,7 +4,7 @@
 
 from Solverz.variable.variables import Vars
 from Solverz.code_printer.python.utilities import *
-
+from Solverz.num_api.module_parser import modules
 
 # %%
 
@@ -162,12 +162,11 @@ def made_numerical(eqs: SymEquations,
                              sparse)
         code['HVP'] = code_HVP
     custom_func = dict()
-    custom_func.update(numerical_interface)
     custom_func.update(parse_trigger_func(eqs.PARAM))
-    F = Solverzlambdify(code_F, 'F_', modules=[custom_func, 'numpy'])
-    J = Solverzlambdify(code_J, 'J_', modules=[custom_func, 'numpy'])
+    F = Solverzlambdify(code_F, 'F_', modules=[custom_func]+modules)
+    J = Solverzlambdify(code_J, 'J_', modules=[custom_func]+modules)
     if make_hvp:
-        HVP = Solverzlambdify(code_HVP, 'Hvp_', modules=[custom_func, 'numpy'])
+        HVP = Solverzlambdify(code_HVP, 'Hvp_', modules=[custom_func]+modules)
     p = parse_p(eqs.PARAM)
     print('Complete!')
     if isinstance(eqs, SymAE) and not isinstance(eqs, SymFDAE):
diff --git a/Solverz/code_printer/python/inline/tests/test_inline_printer.py b/Solverz/code_printer/python/inline/tests/test_inline_printer.py
index b01025d..89bfaa4 100644
--- a/Solverz/code_printer/python/inline/tests/test_inline_printer.py
+++ b/Solverz/code_printer/python/inline/tests/test_inline_printer.py
@@ -5,22 +5,22 @@
 from numpy.testing import assert_allclose
 import re
 import pytest
-from sympy import sin
 from sympy.codegen.ast import Assignment, AddAugmentedAssignment
 
-from Solverz import Model, Var, AliasVar
+from Solverz import Model, Var, AliasVar, sin
 from Solverz.equation.eqn import Eqn, Ode
 from Solverz.equation.equations import DAE, AE
 from Solverz.equation.param import Param
 from Solverz.sym_algebra.symbols import iVar, idx, Para
 from Solverz.variable.variables import combine_Vars, as_Vars
-from Solverz.equation.jac import JacBlock, Ones, Jac
+from Solverz.equation.jac import JacBlock, Jac
 from Solverz.equation.hvp import Hvp
-from Solverz.sym_algebra.functions import Diag, sin, cos, exp
+from Solverz.sym_algebra.functions import *
 from Solverz.code_printer.python.inline.inline_printer import print_J_block, extend, SolList, print_J_blocks, print_J, \
     print_F, print_Hvp, made_numerical, Solverzlambdify
 from Solverz.utilities.address import Address
-from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.module_parser import modules
+
 
 # %%
 row = iVar('row', internal_use=True)
@@ -104,8 +104,8 @@ def test_jbs_printer():
     h = y_[0:1]
     v = y_[1:2]
     g = p_["g"]
-    J_ = zeros((2, 2))
-    J_[0:1,1] += ones(1)
+    J_ = np.zeros((2, 2))
+    J_[0:1,1] += np.ones(1)
     return J_
 """.strip()
 
@@ -118,15 +118,15 @@ def test_jbs_printer():
     data = []
     row.extend([0])
     col.extend([1])
-    data.extend(ones(1))
-    return coo_array((data, (row, col)), (2, 2)).tocsc()
+    data.extend(np.ones(1))
+    return sps.coo_array((data, (row, col)), (2, 2)).tocsc()
 """.strip()
 
 expected_F = """def F_(t, y_, p_):
     h = y_[0:1]
     v = y_[1:2]
     g = p_["g"]
-    _F_ = zeros((2, ))
+    _F_ = np.zeros((2, ))
     _F_[0:1] = v
     _F_[1:2] = g
     return _F_
@@ -141,6 +141,7 @@ def test_print_F_J():
     m.f1 = Ode('f1', f=m.v, diff_var=m.h)
     m.f2 = Ode('f2', f=m.g, diff_var=m.v)
     bball, y0 = m.create_instance()
+    bball.FormJac(y0)
     assert print_J(bball.__class__.__name__,
                    bball.jac,
                    bball.a,
@@ -175,11 +176,11 @@ def test_print_F_J():
     q = y_[3:6]
     p_tag_0 = y_0[0:3]
     q_tag_0 = y_0[3:6]
-    J_ = zeros((6, 6))
-    J_[0:2,1:3] += diagflat(ones(2))
-    J_[2:4,0:2] += diagflat(-ones(2))
-    J_[4:5,2] += -ones(1)
-    J_[5:6,2] += ones(1)
+    J_ = np.zeros((6, 6))
+    J_[0:2,1:3] += np.diagflat(np.ones(2))
+    J_[2:4,0:2] += np.diagflat(-np.ones(2))
+    J_[4:5,2] += -np.ones(1)
+    J_[5:6,2] += np.ones(1)
     return J_
 """.strip()
 
@@ -193,17 +194,17 @@ def test_print_F_J():
     data = []
     row.extend([0, 1])
     col.extend([1, 2])
-    data.extend(ones(2))
+    data.extend(np.ones(2))
     row.extend([2, 3])
     col.extend([0, 1])
-    data.extend(-ones(2))
+    data.extend(-np.ones(2))
     row.extend([4])
     col.extend([2])
-    data.extend(-ones(1))
+    data.extend(-np.ones(1))
     row.extend([5])
     col.extend([2])
-    data.extend(ones(1))
-    return coo_array((data, (row, col)), (6, 6)).tocsc()
+    data.extend(np.ones(1))
+    return sps.coo_array((data, (row, col)), (6, 6)).tocsc()
 """.strip()
 
 expected_F_fdae = """def F_(t, y_, p_, y_0):
@@ -211,7 +212,7 @@ def test_print_F_J():
     q = y_[3:6]
     p_tag_0 = y_0[0:3]
     q_tag_0 = y_0[3:6]
-    _F_ = zeros((6, ))
+    _F_ = np.zeros((6, ))
     _F_[0:2] = -p_tag_0[0:2] + p[1:3]
     _F_[2:4] = p_tag_0[1:3] - p[0:2]
     _F_[4:5] = p_tag_0[0] - p[2]
@@ -237,6 +238,7 @@ def test_print_F_J_FDAE():
     m.ae4 = Eqn('cha4',
                 m.p[2] - m.p0[0])
     fdae, y0 = m.create_instance()
+    fdae.FormJac(y0)
     assert print_J(fdae.__class__.__name__,
                    fdae.jac,
                    fdae.a,
@@ -334,7 +336,7 @@ def test_hvp_printer():
                      var_addr,
                      dict())
 
-    HVP = Solverzlambdify(code, 'Hvp_', [numerical_interface, 'numpy'])
+    HVP = Solverzlambdify(code, 'Hvp_', modules)
     np.testing.assert_allclose(HVP(np.array([1, 2]), dict(), np.array([1, 1])).toarray(),
                                np.array([[2.71828183, -0.90929743],
                                          [0.,  2.]]))
diff --git a/Solverz/code_printer/python/module/module_generator.py b/Solverz/code_printer/python/module/module_generator.py
index 240de04..120eebd 100644
--- a/Solverz/code_printer/python/module/module_generator.py
+++ b/Solverz/code_printer/python/module/module_generator.py
@@ -11,7 +11,7 @@
 from Solverz.code_printer.python.module.module_printer import print_F, print_inner_F, print_sub_inner_F, \
     print_J, print_inner_J, print_Hvp, print_inner_Hvp
 from Solverz.equation.equations import Equations as SymEquations
-from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.module_parser import modules
 from Solverz.variable.variables import Vars, combine_Vars
 from Solverz.equation.hvp import Hvp
 
@@ -67,8 +67,7 @@ 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']
 
-    custom_func = dict()
-    custom_func.update(numerical_interface)
+
 
     def print_trigger_func_code():
         code_tfuc = dict()
@@ -117,7 +116,7 @@ def print_trigger_func_code():
                       p,
                       eqn_parameter,
                       y,
-                      [custom_func, 'numpy'],
+                      modules,
                       numba,
                       directory)
     print('Complete!')
@@ -265,16 +264,24 @@ 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 += 'from numpy import abs\n'
-    code += 'from Solverz.num_api.custom_function import *\n'  # import Solverz built-in func
-    code += 'from scipy.sparse 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 += 'setting = auxiliary["eqn_param"]\n'
     code += 'row = setting["row"]\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 ed965d6..0731404 100644
--- a/Solverz/code_printer/python/module/test/test_module_generator.py
+++ b/Solverz/code_printer/python/module/test/test_module_generator.py
@@ -18,9 +18,15 @@
 from Solverz import load
 auxiliary = load(f"{current_module_dir}\\param_and_setting.pkl")
 from numpy import *
-from numpy import abs
-from Solverz.num_api.custom_function import *
-from scipy.sparse 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
 setting = auxiliary["eqn_param"]
 row = setting["row"]
@@ -99,7 +105,7 @@ def test_AE_module_printer():
         inner_F) == '@njit(cache=True)\ndef inner_F(_F_, x):\n    _F_[0:1] = inner_F0(x)\n    _F_[1:2] = inner_F1(x)\n    return _F_\n'
     assert inspect.getsource(inner_F0.func_code) == '@njit(cache=True)\ndef inner_F0(x):\n    return 2*x[0] + x[1]\n'
     assert inspect.getsource(
-        inner_F1.func_code) == '@njit(cache=True)\ndef inner_F1(x):\n    return x[0]**2 + sin(x[1])\n'
+        inner_F1.func_code) == '@njit(cache=True)\ndef inner_F1(x):\n    return x[0]**2 + np.sin(x[1])\n'
     assert inspect.getsource(
         inner_J) == '@njit(cache=True)\ndef inner_J(_data_, x):\n    _data_[2:3] = inner_J0(x)\n    _data_[3:4] = inner_J1(x)\n    return _data_\n'
 
@@ -125,7 +131,7 @@ def test_AE_module_printer():
     assert inspect.getsource(
         inner_F) == 'def inner_F(_F_, x):\n    _F_[0:1] = inner_F0(x)\n    _F_[1:2] = inner_F1(x)\n    return _F_\n'
     assert inspect.getsource(inner_F0) == 'def inner_F0(x):\n    return 2*x[0] + x[1]\n'
-    assert inspect.getsource(inner_F1) == 'def inner_F1(x):\n    return x[0]**2 + sin(x[1])\n'
+    assert inspect.getsource(inner_F1) == 'def inner_F1(x):\n    return x[0]**2 + np.sin(x[1])\n'
     assert inspect.getsource(
         inner_J) == 'def inner_J(_data_, x):\n    _data_[2:3] = inner_J0(x)\n    _data_[3:4] = inner_J1(x)\n    return _data_\n'
 
@@ -136,7 +142,7 @@ def test_AE_module_printer():
     x = y_[0:2]
     c = p_["c"]
     data_hvp = inner_Hvp(_data_hvp, v_, x, c)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (2, 2)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (2, 2)).tocsc()
 """
 
 expected_inner_Hvp = """@njit(cache=True)
@@ -149,17 +155,17 @@ def inner_Hvp(_data_hvp, v_, x, c):
 
 expected_inner_Hvp0 = """@njit(cache=True)
 def inner_Hvp0(v_, x):
-    return v_[0]*exp(x[0])*ones(1)
+    return v_[0]*np.exp(x[0])*np.ones(1)
 """
 
 expected_inner_Hvp1 = """@njit(cache=True)
 def inner_Hvp1(v_, x):
-    return -v_[1]*sin(x[1])*ones(1)
+    return -v_[1]*sin(x[1])*np.ones(1)
 """
 
 expected_inner_Hvp2 = """@njit(cache=True)
 def inner_Hvp2(c, v_):
-    return v_[1]*c*ones(1)
+    return v_[1]*c*np.ones(1)
 """
 
 
@@ -236,7 +242,7 @@ def inner_F(_F_, p, q, p_tag_0, q_tag_0, pb, qb):
     pb = p_["pb"].get_v_t(t)
     qb = p_["qb"]
     data = inner_J(_data_, p, q, p_tag_0, q_tag_0, pb, qb)
-    return coo_array((data, (row, col)), (164, 164)).tocsc()
+    return sps.coo_array((data, (row, col)), (164, 164)).tocsc()
 """
 
 expected_inner_J = """@njit(cache=True)
@@ -333,7 +339,7 @@ def inner_F(_F_, h, v):
     h = y_[0:1]
     v = y_[1:2]
     data = inner_J(_data_, h, v)
-    return coo_array((data, (row, col)), (2, 2)).tocsc()
+    return sps.coo_array((data, (row, col)), (2, 2)).tocsc()
 """
 
 expected_inner_J1 = """@njit(cache=True)
diff --git a/Solverz/code_printer/python/module/test/test_module_printer.py b/Solverz/code_printer/python/module/test/test_module_printer.py
index ac639ee..0a5a1a9 100644
--- a/Solverz/code_printer/python/module/test/test_module_printer.py
+++ b/Solverz/code_printer/python/module/test/test_module_printer.py
@@ -26,7 +26,7 @@
     lam = p_["lam"]
     ax = ax_trigger_func(x)
     data_hvp = inner_Hvp(_data_hvp, v_, omega, delta, x, y, ax, lam)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
 """.strip()
 expected_FDAE_Hvp = """def Hvp_(t, y_, p_, v_, y_0):
     omega = y_[0:10]
@@ -41,7 +41,7 @@
     lam = p_["lam"]
     ax = ax_trigger_func(x)
     data_hvp = inner_Hvp(_data_hvp, v_, omega, delta, x, y, omega_tag_0, delta_tag_0, x_tag_0, y_tag_0, ax, lam)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
 """.strip()
 
 
@@ -87,13 +87,13 @@ def test_print_Hvp():
     return _data_hvp
 """.strip()
 expected_inner_Hvp0 = """def inner_Hvp0(v_, x):
-    return v_[0]*exp(x[0])*ones(1)
+    return v_[0]*np.exp(x[0])*np.ones(1)
 """.strip()
 expected_inner_Hvp1 = """def inner_Hvp1(v_, x):
-    return -v_[1]*sin(x[1])*ones(1)
+    return -v_[1]*sin(x[1])*np.ones(1)
 """.strip()
 expected_inner_Hvp2 = """def inner_Hvp2(c, v_):
-    return v_[1]*c*ones(1)
+    return v_[1]*c*np.ones(1)
 """.strip()
 
 
@@ -169,7 +169,7 @@ def test_print_inner_Hvp():
     G6 = p_["G6"].get_v_t(t)
     ax = ax_trigger_func(x)
     data = inner_J(_data_, omega, delta, x, y, ax, lam, G6)
-    return coo_array((data, (row, col)), (25, 25)).tocsc()
+    return sps.coo_array((data, (row, col)), (25, 25)).tocsc()
 """.strip()
 expected1 = """def J_(t, y_, p_, y_0):
     omega = y_[0:10]
@@ -185,7 +185,7 @@ def test_print_inner_Hvp():
     G6 = p_["G6"].get_v_t(t)
     ax = ax_trigger_func(x)
     data = inner_J(_data_, omega, delta, x, y, omega_tag_0, delta_tag_0, x_tag_0, y_tag_0, ax, lam, G6)
-    return coo_array((data, (row, col)), (25, 25)).tocsc()
+    return sps.coo_array((data, (row, col)), (25, 25)).tocsc()
 """.strip()
 
 
@@ -234,7 +234,7 @@ def test_print_J():
     return _data_
 """.strip()
 expected3 = """def inner_J0(y):
-    return y*ones(3)
+    return y*np.ones(3)
 """.strip()
 expected4 = """def inner_J1(y):
     return y**2
diff --git a/Solverz/code_printer/python/utilities.py b/Solverz/code_printer/python/utilities.py
index 8706210..b0c386f 100644
--- a/Solverz/code_printer/python/utilities.py
+++ b/Solverz/code_printer/python/utilities.py
@@ -18,11 +18,10 @@
 from Solverz.equation.jac import JacBlock, Jac
 from Solverz.equation.hvp import Hvp
 from Solverz.equation.param import TimeSeriesParam, ParamBase
-from Solverz.sym_algebra.symbols import iVar, SolDict, Para, idx, IdxSymBasic
-from Solverz.sym_algebra.functions import Arange
+from Solverz.sym_algebra.symbols import *
+from Solverz.sym_algebra.functions import *
 from Solverz.utilities.address import Address
 from Solverz.utilities.type_checker import is_number
-from Solverz.num_api.custom_function import numerical_interface
 from Solverz.num_api.num_eqn import nAE, nFDAE, nDAE
 
 
@@ -193,71 +192,5 @@ def print_eqn_assignment(EQNs: Dict[str, Eqn],
     return eqn_declaration
 
 
-class coo_2_csc(Symbol):
 
-    def __new__(cls, eqn_size: int, vsize: int):
-        obj = Symbol.__new__(cls, f'coo_2_csc')
-        obj.eqn_size = eqn_size
-        obj.vsize = vsize
-        return obj
 
-    def _numpycode(self, printer, **kwargs):
-        return f'coo_array((data, (row, col)), ({self.eqn_size}, {self.vsize})).tocsc()'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-
-
-class coo_2_csc_hvp(Symbol):
-
-    def __new__(cls, eqn_size: int, vsize: int):
-        obj = Symbol.__new__(cls, f'coo_2_csc')
-        obj.eqn_size = eqn_size
-        obj.vsize = vsize
-        return obj
-
-    def _numpycode(self, printer, **kwargs):
-        return f'coo_array((data_hvp, (row_hvp, col_hvp)), ({self.eqn_size}, {self.vsize})).tocsc()'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-
-
-class coo_array(Function):
-
-    @classmethod
-    def eval(cls, *args):
-        if len(args) > 1:
-            raise ValueError(
-                f"Solverz' coo_array object accepts only one inputs.")
-
-    def _numpycode(self, printer, **kwargs):
-        return f'coo_array({printer._print(self.args[0])})'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-
-
-class extend(Function):
-
-    def _numpycode(self, printer, **kwargs):
-        return f'{printer._print(self.args[0])}.extend({printer._print(self.args[1])})'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-
-
-class zeros(Function):
-    # print zeros(6,6) as zeros((6,6))
-    # or zeros(6,) as zeros((6,))
-    def _numpycode(self, printer, **kwargs):
-        if len(self.args) == 2:
-            temp1 = printer._print(self.args[0])
-            temp2 = printer._print(self.args[1])
-            return r'zeros((' + temp1 + ', ' + temp2 + r'))'
-        elif len(self.args) == 1:
-            temp = printer._print(self.args[0])
-            return r'zeros((' + temp + ', ' + r'))'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
diff --git a/Solverz/equation/eqn.py b/Solverz/equation/eqn.py
index 640b99c..1e7c061 100644
--- a/Solverz/equation/eqn.py
+++ b/Solverz/equation/eqn.py
@@ -3,27 +3,19 @@
 from typing import Union, List, Dict, Callable
 
 import numpy as np
-from sympy import Symbol, Expr, latex, Derivative, sympify
+from sympy import Symbol, Expr, latex, Derivative, sympify, Function
 from sympy import lambdify as splambdify
 from sympy.abc import t, x
 
 from Solverz.sym_algebra.symbols import iVar, Para, IdxVar, idx, IdxPara, iAliasVar, IdxAliasVar
 from Solverz.variable.ssymbol import Var
-from Solverz.sym_algebra.functions import Mat_Mul, Slice, F
+from Solverz.sym_algebra.functions import Mat_Mul, Slice
 from Solverz.sym_algebra.matrix_calculus import MixedEquationDiff
-from Solverz.sym_algebra.transform import finite_difference, semi_descritize
-from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.module_parser import modules
 from Solverz.variable.ssymbol import sSym2Sym
 from Solverz.utilities.type_checker import is_zero
 
 
-def sVar2Var(var: Union[Var, iVar, List[iVar, Var]]) -> Union[iVar, List[iVar]]:
-    if isinstance(var, list):
-        return [arg.symbol if isinstance(arg, Var) else arg for arg in var]
-    else:
-        return var.symbol if isinstance(var, Var) else var
-
-
 class Eqn:
     """
     The Equation object
@@ -36,7 +28,7 @@ def __init__(self,
             raise ValueError("Equation name must be string!")
         self.name: str = name
         self.LHS = 0
-        self.RHS = sympify(eqn)
+        self.RHS = sympify(sSym2Sym(eqn))
         self.SYMBOLS: Dict[str, Symbol] = self.obtain_symbols()
 
         # if the eqn has Mat_Mul, then label it as mixed-matrix-vector equation
@@ -62,7 +54,7 @@ def obtain_symbols(self) -> Dict[str, Symbol]:
         return sorted_dict
 
     def lambdify(self) -> Callable:
-        return splambdify(self.SYMBOLS.values(), self.RHS, [numerical_interface, 'numpy'])
+        return splambdify(self.SYMBOLS.values(), self.RHS, modules)
 
     def eval(self, *args: Union[np.ndarray]) -> np.ndarray:
         return self.NUM_EQN(*args)
@@ -134,7 +126,7 @@ def __init__(self, name: str, eqn: Expr, diff_var: Symbol, var_idx=None):
                 self.var_idx_func = Eqn('To evaluate var_idx of variable' + self.diff_var.name, Slice(*temp))
             elif isinstance(self.var_idx, Expr):
                 self.var_idx_func = Eqn('To evaluate var_idx of variable' + self.diff_var.name, self.var_idx)
-        self.LHS = Derivative(F, diff_var)
+        self.LHS = Derivative(Function('F'), diff_var)
         self.dim = -1
         self.v_type = ''
 
@@ -154,340 +146,7 @@ class Ode(Eqn):
     def __init__(self, name: str,
                  f,
                  diff_var: Union[iVar, IdxVar, Var]):
-        super().__init__(name, sSym2Sym(f))
+        super().__init__(name, f)
         diff_var = sSym2Sym(diff_var)
         self.diff_var = diff_var
         self.LHS = Derivative(diff_var, t)
-
-
-class Pde(Eqn):
-    """
-    The class of partial differential equations
-    """
-    pass
-
-
-class HyperbolicPde(Pde):
-    r"""
-    The class for hyperbolic PDE reading
-
-    .. math::
-
-         \frac{\partial{u}}{\partial{t}}+\frac{\partial{f(u)}}{\partial{x}}=S(u)
-
-    where $u$ is the state vector, $f(u)$ is the flux function and $S(u)$ is the source term.
-
-    Parameters
-    ==========
-
-    two_dim_var : iVar or list of Var
-
-        Specify the two-dimensional variables in the PDE. Some of the variables, for example, the mass flow $\dot{m}$ in
-        the heat transmission equation, are not two-dimensional variables.
-
-    """
-
-    def __init__(self, name: str,
-                 diff_var: iVar | Var,
-                 flux: Expr = 0,
-                 source: Expr = 0,
-                 two_dim_var: Union[iVar, Var, List[iVar | Var]] = None):
-        if isinstance(source, (float, int)):
-            source = sympify(source)
-        super().__init__(name, source)
-        diff_var = sVar2Var(diff_var)
-        two_dim_var = sVar2Var(two_dim_var) if two_dim_var is not None else None
-        self.diff_var = diff_var
-        if isinstance(flux, (float, int)):
-            flux = sympify(flux)
-        if isinstance(source, (float, int)):
-            flux = sympify(flux)
-        self.flux = flux
-        self.source = source
-        self.two_dim_var = [two_dim_var] if isinstance(two_dim_var, iVar) else two_dim_var
-        self.LHS = Derivative(diff_var, t) + Derivative(flux, x)
-
-    def derive_derivative(self):
-        pass
-
-    def finite_difference(self, scheme='central diff', direction=None, M: int = 0, dx=None):
-        r"""
-        Discretize hyperbolic PDE as AEs.
-
-        Parameters
-        ==========
-
-        scheme : str
-
-            1 - Central difference
-
-            .. math::
-
-                \frac{\partial{u}}{\partial{t}}\approx\frac{u_{i+1}^{j+1}-u_{i+1}^{j}+u_{i}^{j+1}-u_{i}^{j}}{2\Delta t}
-
-            .. math::
-
-                \frac{\partial{f(u)}}{\partial{x}}\approx\frac{f(u_{i+1}^{j+1})-f(u_{i}^{j+1})+f(u_{i+1}^{j})-f(u_{i}^{j})}{2\Delta x}
-
-            .. math::
-
-                S(u)\approx S\left(\frac{u_{i+1}^{j+1}+u_{i}^{j+1}+u_{i+1}^{j}+u_{i}^{j}}{4}\right)
-
-            2 - Backward Time Backward/Forward Space
-
-            If direction equals 1, then do backward space difference, which derives
-
-            .. math::
-
-                \frac{\partial{u}}{\partial{t}}\approx\frac{u_{i+1}^{j+1}-u_{i+1}^{j}}{\Delta t}
-
-            .. math::
-
-                \frac{\partial{f(u)}}{\partial{x}}\approx\frac{f(u_{i+1}^{j+1})-f(u_{i}^{j+1})}{\Delta x}
-
-            .. math::
-
-                S(u)\approx S\left(u_{i+1}^{j+1}\right)
-
-            If direction equals -1, then do forward space difference, which derives
-
-            .. math::
-
-                \frac{\partial{u}}{\partial{t}}\approx\frac{u_{i}^{j+1}-u_{i}^{j}}{\Delta t}
-
-            .. math::
-
-                \frac{\partial{f(u)}}{\partial{x}}\approx\frac{f(u_{i+1}^{j+1})-f(u_{i}^{j+1})}{\Delta x}
-
-            .. math::
-
-                S(u)\approx S\left(u_{i}^{j+1}\right)
-
-        direction : int
-
-            To tell which side of boundary conditions is given in scheme 2.
-
-        M : int
-
-        The total number of spatial sections.
-
-        dx : Number
-
-        Spatial difference step size
-
-        Returns
-        =======
-
-        AE : Eqn
-
-            Let's take central difference as an example, this function returns the algebraic equation
-
-            .. math::
-
-                \begin{aligned}
-                    0=&\Delta x(\tilde{u}[1:M]-\tilde{u}^0[1:M]+\tilde{u}[0:M-1]-\tilde{u}^0[0:M-1])+\\
-                      &\Delta t(f(\tilde{u}[1:M])-f(\tilde{u}[0:M-1])+f(\tilde{u}^0[1:M])-f(\tilde{u}^0[0:M-1]))+\\
-                      &2\Delta x\Delta t\cdot S\left(\tilde{u}[1:M]-\tilde{u}^0[1:M]+\tilde{u}[0:M-1]-\tilde{u}^0[0:M-1]}{4}\right)
-                \end{aligned}
-
-            where we denote by vector $\tilde{u}$ the discrete spatial distribution of state $u$, by $\tilde{u}^0$ the
-            initial value of $\tilde{u}$, and by $M$ the last index of $\tilde{u}$.
-
-        """
-        if isinstance(M, (int, np.integer)):
-            if M < 0:
-                raise ValueError(f'Total nunmber of PDE sections {M} < 0')
-        else:
-            raise TypeError(f'Do not support M of type {type(M)}')
-        if M == 0:
-            M = idx('M')
-
-        if scheme == 'central diff':
-            return Eqn('FDM of ' + self.name + 'w.r.t.' + self.diff_var.name + 'using central diff',
-                       finite_difference(self.diff_var,
-                                         self.flux,
-                                         self.source,
-                                         self.two_dim_var,
-                                         M,
-                                         'central diff',
-                                         dx=dx))
-        elif scheme == 'euler':
-            return Eqn('FDM of ' + self.name + 'w.r.t.' + self.diff_var.name + 'using Euler',
-                       finite_difference(self.diff_var,
-                                         self.flux,
-                                         self.source,
-                                         self.two_dim_var,
-                                         M,
-                                         'euler',
-                                         direction=direction,
-                                         dx=dx))
-
-    def semi_discretize(self,
-                        a0=None,
-                        a1=None,
-                        scheme='TVD1',
-                        M: int = 0,
-                        output_boundary=True,
-                        dx=None) -> List[Eqn]:
-        r"""
-        Semi-discretize the hyperbolic PDE of nonlinear conservation law as ODEs using the Kurganov-Tadmor scheme
-        (see [Kurganov2000]_). The difference stencil is as follows, with $x_{j+1}-x_{j}=\Delta x$.
-
-            .. image:: ../../pics/difference_stencil.png
-               :height: 100
-
-        Parameters
-        ==========
-
-        a0 : Expr
-
-            Maximum local speed $a_{j+1/2}$, with formula
-
-            .. math::
-
-                a_{j+1/2}=\max\qty{\rho\qty(\pdv{f}{u}\qty(u^+_{j+1/2})),\rho\qty(\pdv{f}{u}\qty(u^-_{j+1/2}))},
-
-            where
-
-            .. math::
-
-                \rho(A)=\max_i|\lambda_i(A)|.
-
-            If $a_0$ or $a_1$ is None, then they will be set as ``Para`` ``ajp12`` and ``ajm12`` respectively.  
-
-        a1 : Expr
-
-            Maximum local speed $a_{j-1/2}$, with formula
-
-            .. math::
-
-                a_{j-1/2}=\max\qty{\rho\qty(\pdv{f}{u}\qty(u^+_{j-1/2})),\rho\qty(\pdv{f}{u}\qty(u^-_{j-1/2}))}.
-
-        scheme : str
-
-            If scheme==1, 2nd scheme else, else, use 1st scheme.
-
-        M : int
-
-            The total number of spatial sections.
-
-        output_boundary : bool
-
-            If true, output equations about the boundary conditions. For example,
-
-           >>> from Solverz import HyperbolicPde, iVar
-           >>> T = iVar('T')
-           >>> p = HyperbolicPde(name = 'heat transfer', diff_var=T, flux=T)
-           >>> p.semi_discretize(a0=1,a2=1, scheme=2, M=2, output_boundary=True)
-           1
-           >>> p.semi_discretize(a0=1,a2=1, scheme=2, M=2, output_boundary=False)
-           2
-
-        dx : Number
-
-            spatial difference step size
-
-        Returns
-        =======
-
-        ODE : List[Union[Ode, Eqn]]
-
-            This function returns the for $2\leq j\leq M-2$
-
-            .. math::
-
-                \dv{t}u_j=-\frac{H_{j+1/2}-H_{j-1/2}}{\Delta x}+S(u_j)
-
-            and for $j=1,M-1$
-
-            .. math::
-
-                \dv{t}u_j=-\frac{f(u_{j+1})-f(u_{j-1})}{2\Delta x}+\frac{a_{j+1/2}(u_{j+1}-u_j)-a_{j-1/2}(u_j-u_{j-1})}{2\Delta x}+S(u_j),
-
-            where
-
-            .. math::
-
-                H_{j+1/2}=\frac{f(u^+_{j+1/2})+f(u^-_{j+1/2})}{2}-\frac{a_{j+1/2}}{2}\qty[u^+_{j+1/2}-u^-_{j+1/2}],
-
-            .. math::
-
-                H_{j-1/2}=\frac{f(u^+_{j-1/2})+f(u^-_{j-1/2})}{2}-\frac{a_{j-1/2}}{2}\qty[u^+_{j-1/2}-u^-_{j-1/2}],
-
-            .. math::
-
-                u^+_{j+1/2}=u_{j+1}-\frac{\Delta x}{2}(u_x)_{j+1},\quad u^-_{j+1/2}=u_j+\frac{\Delta x}{2}(u_x)_j,
-
-            .. math::
-
-                u^+_{j-1/2}=u_{j}-\frac{\Delta x}{2}(u_x)_{j},\quad u^-_{j-1/2}=u_{j-1}+\frac{\Delta x}{2}(u_x)_{j-1},
-
-            .. math::
-
-                (u_x)_j=\operatorname{minmod}\qty(\theta\frac{u_j-u_{j-1}}{\Delta x},\frac{u_{j+1}-u_{j-1}}{2\Delta x},\theta\frac{u_{j+1}-u_{j}}{\Delta x}),\quad \theta\in[1,2],
-
-            and by linear extrapolation
-
-            .. math::
-
-                u_0=2u_\text{L}-u_1,\quad u_M=2u_\text{R}-u_{M-1}.
-
-
-        .. [Kurganov2000] Alexander Kurganov, Eitan Tadmor, New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations, Journal of Computational Physics, Volume 160, Issue 1, 2000, Pages 241-282, `<https://doi.org/10.1006/jcph.2000.6459>`_
-
-        """
-        if isinstance(M, (int, np.integer)):
-            if M < 0:
-                raise ValueError(f'Total nunmber of PDE sections {M} < 0')
-        else:
-            raise TypeError(f'Do not support M of type {type(M)}')
-        if M == 0:
-            M = idx('M')
-        u = self.diff_var
-        dae_list = []
-        if scheme == 'TVD2':
-            eqn_dict = semi_descritize(self.diff_var,
-                                       self.flux,
-                                       self.source,
-                                       self.two_dim_var,
-                                       M,
-                                       scheme='TVD2',
-                                       a0=a0,
-                                       a1=a1,
-                                       dx=dx)
-            dae_list.extend([Ode('SDM of ' + self.name + ' 1',
-                                 eqn_dict['Ode'][0][0],
-                                 eqn_dict['Ode'][0][1]),
-                             Ode('SDM of ' + self.name + ' 2',
-                                 eqn_dict['Ode'][1][0],
-                                 eqn_dict['Ode'][1][1]),
-                             Ode('SDM of ' + self.name + ' 3',
-                                 eqn_dict['Ode'][2][0],
-                                 eqn_dict['Ode'][2][1]),
-                             Eqn('minmod limiter 1 of ' + u.name,
-                                 eqn_dict['Eqn'][0]),
-                             Eqn('minmod limiter 2 of ' + u.name,
-                                 eqn_dict['Eqn'][1]),
-                             Eqn('minmod limiter 3 of ' + u.name,
-                                 eqn_dict['Eqn'][2])])
-        elif scheme == 'TVD1':
-            eqn_dict = semi_descritize(self.diff_var,
-                                       self.flux,
-                                       self.source,
-                                       self.two_dim_var,
-                                       M,
-                                       scheme='TVD1',
-                                       a0=a0,
-                                       a1=a1,
-                                       dx=dx)
-            dae_list.extend([Ode('SDM of ' + self.name + 'using',
-                                 eqn_dict['Ode'][0],
-                                 eqn_dict['Ode'][1])])
-        else:
-            raise NotImplementedError(f'Scheme {scheme} not implemented')
-
-        if output_boundary:
-            dae_list.extend([Eqn(f'Equation of {u[M]}', u[M] - 2 * iVar(u.name + 'R') + u[M - 1]),
-                             Eqn(f'Equation of {u[0]}', u[0] - 2 * iVar(u.name + 'L') + u[1])])
-
-        return dae_list
diff --git a/Solverz/equation/jac.py b/Solverz/equation/jac.py
index 5ff84ba..ced90ef 100644
--- a/Solverz/equation/jac.py
+++ b/Solverz/equation/jac.py
@@ -6,7 +6,8 @@
 from sympy import Expr, Function, Integer
 from Solverz.sym_algebra.symbols import iVar, IdxVar
 from Solverz.utilities.type_checker import is_vector, is_scalar, is_integer, is_number, PyNumber, is_zero
-from Solverz.sym_algebra.functions import Diag
+from Solverz.sym_algebra.functions import Diag, Ones
+
 
 SolVar = Union[iVar, IdxVar]
 
@@ -297,27 +298,5 @@ def slice2array(s: slice) -> np.ndarray:
     return np.arange(s.start, s.stop)
 
 
-class Ones(Function):
-    r"""
-    The all-one vector to broadcast scalars to a vector
-    For example, 2 -> 2*Ones(eqn_size)
-    The derivative is d(Ones(eqn_size))/dx=0
-    """
-
-    @classmethod
-    def eval(cls, x):
-        if not isinstance(x, Integer):
-            raise ValueError("The arg of Ones() should be integer.")
-
-    def _eval_derivative(self, s):
-        return 0
-
-    def _numpycode(self, printer, **kwargs):
-        x = self.args[0]
-        return r'ones(' + printer._print(x, **kwargs) + r')'
 
-    def _lambdacode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
 
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
diff --git a/Solverz/equation/test/test_eqn.py b/Solverz/equation/test/test_eqn.py
index bb21504..630092c 100644
--- a/Solverz/equation/test/test_eqn.py
+++ b/Solverz/equation/test/test_eqn.py
@@ -22,11 +22,24 @@
 assert isinstance(f1, np.ndarray)
 assert f1.ndim == 1
 
-# discard zero derivative
-u = Var('u', [1, 1, 1])
-e = Var('e', [1, 1, 1])
-umin = Param('umin', [0,0,0])
-umax = Param('umax', [5,5,5])
-F = Ode('Anti', AntiWindUp(u, umin, umax, e), u)
-F.derive_derivative()
-assert 'u' not in F.derivatives
+
+def test_discarding_zero_deri():
+    # discard zero derivative
+    u = Var('u', [1, 1, 1])
+    e = Var('e', [1, 1, 1])
+    umin = Param('umin', [0,0,0])
+    umax = Param('umax', [5,5,5])
+    F = Ode('Anti', AntiWindUp(u, umin, umax, e), u)
+    F.derive_derivative()
+    assert 'u' not in F.derivatives
+
+def test_Var_converter():
+    # convert Var/Param to iVar/Para objects
+    f = Eqn('f', Var('x'))
+    assert f.SYMBOLS['x'].__str__() == 'x'
+    f = Eqn('f', Param('x'))
+    assert f.SYMBOLS['x'].__str__() == 'x'
+    g = Ode('g', Var('x'), diff_var=Var('y'))
+    assert g.SYMBOLS['x'].__str__() == 'x'
+    assert g.diff_var.name == 'y'
+
diff --git a/Solverz/model/basic.py b/Solverz/model/basic.py
index ea0ff2c..eaffb0e 100644
--- a/Solverz/model/basic.py
+++ b/Solverz/model/basic.py
@@ -85,7 +85,8 @@ def create_instance(self):
         if eqn_type == 'FDAE':
             for i in range(nstep):
                 eqs.update_param(y0.derive_alias(f'_tag_{i}'))
-        eqs.FormJac(y0)
+        # eqs.FormJac(y0)
+        eqs.assign_eqn_var_address(y0)
 
         if eqs.eqn_size != eqs.vsize:
             warnings.warn(f'Equation size {eqs.eqn_size} and variable size {eqs.vsize} not equal!')
diff --git a/Solverz/num_api/custom_function.py b/Solverz/num_api/custom_function.py
index 6cbd9dc..2ffa525 100644
--- a/Solverz/num_api/custom_function.py
+++ b/Solverz/num_api/custom_function.py
@@ -2,28 +2,17 @@
 
 from functools import reduce
 
+import warnings
 import numpy as np
 from numpy import linalg
 from scipy.sparse import diags, csc_array, coo_array, linalg as sla
 from numba import njit
 
+
 # from cvxopt.umfpack import linsolve
 # from cvxopt import matrix, spmatrix
 
-numerical_interface = {'coo_array': coo_array, 'csc_array': csc_array}
-
-
-def implements_nfunc(nfunc_name: str):
-    """Register an DT function implementation for sympy Expr."""
-
-    def decorator(func):
-        numerical_interface[nfunc_name] = func
-        return func
-
-    return decorator
-
 
-@implements_nfunc('sol_slice')
 def sol_slice(*args):
     """
     This is used to convert the slice arguments to int
@@ -31,7 +20,6 @@ def sol_slice(*args):
     return slice(*[int(arg_[0]) if isinstance(arg_, np.ndarray) else arg_ for arg_ in args])
 
 
-@implements_nfunc('Slice')
 def Slice(*args):
     """
     This is used to evaluate the slice index of IdxVar/IdxParam/IdxConst
@@ -39,57 +27,20 @@ def Slice(*args):
     return sol_slice(*args)
 
 
-@implements_nfunc('ix_')
 def ix_(arg: np.ndarray):
     return arg.reshape(-1, )
 
 
-@implements_nfunc('Sign')
 def _sign(arg):
     return np.sign(arg)
 
 
-@implements_nfunc('minmod')
-def minmod(a, b, c):
-    stacked_array = np.hstack((a, b, c))
-    cd1 = (a > 0) & (b > 0) & (c > 0)
-    cd2 = (a < 0) & (b < 0) & (c < 0)
-    conditions = [cd1,
-                  cd2]
-    choice_list = [np.min(stacked_array, axis=0),
-                   np.max(stacked_array, axis=0)]
-    return np.select(conditions, choice_list, 0)
-
-
-@implements_nfunc('minmod_flag')
-def minmod_flag(*args):
-    if len(args) != 3:
-        raise ValueError("Input arg length must be 3")
-
-    shapes = [arg.shape[0] for arg in args]
-
-    a, b, c = args
-    stacked_array = np.hstack((a, b, c))
-    if all(x == shapes[0] for x in shapes):
-        cd1 = (a > 0) & (b > 0) & (c > 0)
-        cd2 = (a < 0) & (b < 0) & (c < 0)
-        conditions = [cd1,
-                      cd2]
-        choice_list = [np.argmin(stacked_array, axis=0),
-                       np.argmax(stacked_array, axis=0)]
-    else:
-        raise ValueError(f"Length of Input array not consistent {shapes}")
-    return np.select(conditions, choice_list, 3)
-
-
-@implements_nfunc('Heaviside')
-@njit
+@njit(cache=True)
 def Heaviside(x):
     return np.where(x >= 0, 1.0, 0.0)
 
 
-@implements_nfunc('switch')
-@njit
+@njit(cache=True)
 def switch(*args):
     flag = args[-1]
     flag_shape = args[-1].shape
@@ -111,49 +62,48 @@ def switch(*args):
     return np.select(conditions, choice_list, 0)
 
 
-@implements_nfunc('SolIn')
 @njit(cache=True)
-def SolIn(x, xmin, xmax):
+def Saturation(x, xmin, xmax):
+    x = np.asarray(x).reshape((-1,))
+    return np.clip(x, xmin, xmax)
+
+
+@njit(cache=True)
+def In(x, xmin, xmax):
     x = np.asarray(x).reshape((-1,))
     return np.bitwise_and(x >= xmin, x <= xmax).astype(np.int32)
 
 
-@implements_nfunc('SolGreaterThan')
 @njit(cache=True)
-def SolGreaterThan(x, y):
+def GreaterThan(x, y):
     x = np.asarray(x).reshape((-1,))
     return (x > y).astype(np.int32)
 
 
-@implements_nfunc('SolLessThan')
 @njit(cache=True)
-def SolLessThan(x, y):
+def LessThan(x, y):
     x = np.asarray(x).reshape((-1,))
     return (x < y).astype(np.int32)
 
 
-@implements_nfunc('And')
 @njit(cache=True)
 def And(x, y):
     x = np.asarray(x).reshape((-1,))
     return x & y
 
 
-@implements_nfunc('Or')
 @njit(cache=True)
 def Or(x, y):
     x = np.asarray(x).reshape((-1,))
     return x | y
 
 
-@implements_nfunc('Not')
 @njit(cache=True)
 def Not(x):
     x = np.asarray(x).reshape((-1,))
     return np.ones_like(x) - x
 
 
-@implements_nfunc('Diag')
 def diag(x) -> np.ndarray:
     """
     Generate diagonal matrix of given vector X
@@ -166,7 +116,6 @@ def diag(x) -> np.ndarray:
         return np.diagflat(x)
 
 
-@implements_nfunc('dConv_s')
 def DT_conv(*args, method='conv') -> np.ndarray:
     r"""
     Perform the convolutions in DT computations.
@@ -205,7 +154,6 @@ def DT_conv(*args, method='conv') -> np.ndarray:
         return np.array(np.real(np.fft.ifft(reduce(lambda a, b: a * b, y))[k - 1]))
 
 
-@implements_nfunc('dLinspace')
 def linspace(start, end) -> np.ndarray:
     r"""
 
diff --git a/Solverz/num_api/module_parser.py b/Solverz/num_api/module_parser.py
new file mode 100644
index 0000000..7c0fd44
--- /dev/null
+++ b/Solverz/num_api/module_parser.py
@@ -0,0 +1,16 @@
+# parse modules from built-in custom functions
+import numpy
+import scipy
+import warnings
+
+import Solverz.num_api.custom_function as SolCF
+
+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
+
+# parse modules from museum
+try:
+    import SolMuseum.num_api as SolMF
+    modules[0]['SolMF'] = SolMF
+except ModuleNotFoundError as e:
+    warnings.warn(f'Failed to import num api from SolMuseum: {e}')
diff --git a/Solverz/num_api/test/test_custom_func.py b/Solverz/num_api/test/test_custom_func.py
index d853929..3cef960 100644
--- a/Solverz/num_api/test/test_custom_func.py
+++ b/Solverz/num_api/test/test_custom_func.py
@@ -1,31 +1,31 @@
 import numpy as np
 
-from Solverz.num_api.custom_function import SolIn, SolLessThan, SolGreaterThan, And, Or, Not
+from Solverz.num_api.custom_function import In, LessThan, GreaterThan, And, Or, Not
 
 
 def test_in():
     a = np.array([1, 2, 3.0])
     b = np.array([3, 2, 1.0])
-    np.testing.assert_allclose(SolIn(a, 2, 2), np.array([0, 1, 0]))
-    np.testing.assert_allclose(SolIn(a, 1, 3), np.array([1, 1, 1]))
-    np.testing.assert_allclose(SolIn(a, 1, 2.5), np.array([1, 1, 0]))
-    np.testing.assert_allclose(SolIn(5, 4, 4.5), np.array([0]))
+    np.testing.assert_allclose(In(a, 2, 2), np.array([0, 1, 0]))
+    np.testing.assert_allclose(In(a, 1, 3), np.array([1, 1, 1]))
+    np.testing.assert_allclose(In(a, 1, 2.5), np.array([1, 1, 0]))
+    np.testing.assert_allclose(In(5, 4, 4.5), np.array([0]))
 
 
 def test_gt():
     a = np.array([1, 2, 3.0])
     b = np.array([3, 2, 1.0])
-    np.testing.assert_allclose(SolGreaterThan(a, b), np.array([0, 0, 1]))
-    np.testing.assert_allclose(SolGreaterThan(5, 4), np.array([1]))
-    np.testing.assert_allclose(SolGreaterThan(5, 6), np.array([0]))
+    np.testing.assert_allclose(GreaterThan(a, b), np.array([0, 0, 1]))
+    np.testing.assert_allclose(GreaterThan(5, 4), np.array([1]))
+    np.testing.assert_allclose(GreaterThan(5, 6), np.array([0]))
 
 
 def test_lt():
     a = np.array([1, 2, 3.0])
     b = np.array([3, 2, 1.0])
-    np.testing.assert_allclose(SolLessThan(a, b), np.array([1, 0, 0]))
-    np.testing.assert_allclose(SolLessThan(5, 4), np.array([0]))
-    np.testing.assert_allclose(SolLessThan(5, 6), np.array([1]))
+    np.testing.assert_allclose(LessThan(a, b), np.array([1, 0, 0]))
+    np.testing.assert_allclose(LessThan(5, 4), np.array([0]))
+    np.testing.assert_allclose(LessThan(5, 6), np.array([1]))
 
 
 def test_and():
diff --git a/Solverz/sym_algebra/functions.py b/Solverz/sym_algebra/functions.py
index 3133b6e..5006963 100644
--- a/Solverz/sym_algebra/functions.py
+++ b/Solverz/sym_algebra/functions.py
@@ -6,14 +6,6 @@
 from Solverz.variable.ssymbol import sSym2Sym
 
 
-# %% miscellaneous
-class F(Function):
-    """
-    For the usage of denoting the function being differentiated in EqnDiff object only
-    """
-    pass
-
-
 # %%
 def VarParser(cls):
     # To convert non-sympy symbols to sympy symbols
@@ -159,7 +151,7 @@ def _latex(self, printer, **kwargs):
 
     def _numpycode(self, printer, **kwargs):
 
-        return r'diagflat(' + printer._print(self.args[0], **kwargs) + r')'
+        return r'np.diagflat(' + printer._print(self.args[0], **kwargs) + r')'
 
     def _lambdacode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
@@ -174,10 +166,21 @@ def _octave(self, printer, **kwargs):
 # %% Univariate func
 @VarParser
 class UniVarFunc(Function):
+    arglength = 1
+
     @classmethod
     def eval(cls, *args):
         if len(args) != 1:
-            raise TypeError(f'Supports one operand while {len(args)} input!')
+            raise TypeError(f'{cls.name} supports {cls.arglength} operand while {len(args)} input!')
+
+    def _numpycode(self, printer, **kwargs):
+        pass
+
+    def _lambdacode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
 
 
 class Abs(UniVarFunc):
@@ -206,10 +209,7 @@ def fdiff(self, argindex=1):
             raise ArgumentIndexError(self, argindex)
 
     def _numpycode(self, printer, **kwargs):
-        return r'abs(' + printer._print(self.args[0]) + r')'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.abs(' + printer._print(self.args[0]) + r')'
 
 
 class exp(UniVarFunc):
@@ -221,10 +221,7 @@ def fdiff(self, argindex=1):
         return exp(*self.args)
 
     def _numpycode(self, printer, **kwargs):
-        return r'exp(' + printer._print(self.args[0]) + r')'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.exp(' + printer._print(self.args[0]) + r')'
 
 
 class ln(UniVarFunc):
@@ -236,12 +233,11 @@ def fdiff(self, argindex=1):
         return 1 / self.args[0]
 
     def _numpycode(self, printer, **kwargs):
-        return r'log(' + printer._print(self.args[0]) + r')'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.log(' + printer._print(self.args[0]) + r')'
 
 
+# Notice: Do not succeed sympy.sin or cos here because the args of symbolic functions are Solverz.ssymbol.
+# We have to parse them first.
 class sin(UniVarFunc):
     r"""
     The sine function.
@@ -254,10 +250,7 @@ def fdiff(self, argindex=1):
             raise ArgumentIndexError(self, argindex)
 
     def _numpycode(self, printer, **kwargs):
-        return r'sin(' + printer._print(self.args[0]) + r')'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.sin(' + printer._print(self.args[0]) + r')'
 
 
 class cos(UniVarFunc):
@@ -272,10 +265,7 @@ def fdiff(self, argindex=1):
             raise ArgumentIndexError(self, argindex)
 
     def _numpycode(self, printer, **kwargs):
-        return r'cos(' + printer._print(self.args[0]) + r')'
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.cos(' + printer._print(self.args[0]) + r')'
 
 
 class Sign(UniVarFunc):
@@ -300,16 +290,10 @@ def fdiff(self, argindex=1):
         raise ArgumentIndexError(self, argindex)
 
     def _numpycode(self, printer, **kwargs):
-        return r'sign(' + printer._print(self.args[0], **kwargs) + r')'
-
-    def _lambdacode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
-
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+        return r'np.sign(' + printer._print(self.args[0], **kwargs) + r')'
 
     def _octave(self, printer, **kwargs):
-        return r'sign(' + printer._print(self.args[0], **kwargs) + r')'
+        return r'np.sign(' + printer._print(self.args[0], **kwargs) + r')'
 
 
 class heaviside(UniVarFunc):
@@ -334,41 +318,40 @@ def fdiff(self, argindex=1):
         raise ArgumentIndexError(self, argindex)
 
     def _numpycode(self, printer, **kwargs):
+        return r'SolCF.Heaviside(' + printer._print(self.args[0], **kwargs) + r')'
+
+    def _octave(self, printer, **kwargs):
         return r'Heaviside(' + printer._print(self.args[0], **kwargs) + r')'
 
-    def _lambdacode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
 
-    def _pythoncode(self, printer, **kwargs):
-        return self._numpycode(printer, **kwargs)
+class Not(UniVarFunc):
+    """
+    Represents bitwise_not
+    """
 
-    def _octave(self, printer, **kwargs):
-        return r'Heaviside(' + printer._print(self.args[0], **kwargs) + r')'
+    def _eval_derivative(self, s):
+        return Integer(0)
+
+    def _sympystr(self, printer, **kwargs):
+        return '({op1})'.format(op1=printer._print(self.args[0]))
+
+    def _numpycode(self, printer, **kwargs):
+        return r'SolCF.Not(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 # %% multi-variate func
 @VarParser
 class MulVarFunc(Function):
-    pass
-
-
-class minmod_flag(MulVarFunc):
-    """
-    Different from `minmod`, minmod function outputs the position of args instead of the values of args.
-    """
+    arglength = 3
 
     @classmethod
     def eval(cls, *args):
-        if len(args) != 3:
-            raise TypeError(f"minmod takes 3 positional arguments but {len(args)} were given!")
-
-
-class switch(MulVarFunc):
-    def _eval_derivative(self, s):
-        return switch(*[arg.diff(s) for arg in self.args[0:len(self.args) - 1]], self.args[-1])
+        if len(args) != cls.arglength:
+            raise TypeError(f"{cls.__name__} takes {cls.arglength} positional arguments but {len(args)} were given!")
 
     def _numpycode(self, printer, **kwargs):
-        return r'switch(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
+        return (f'{self.__class__.__name__}' + r'(' +
+                ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')')
 
     def _lambdacode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
@@ -377,6 +360,14 @@ def _pythoncode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
 
 
+class switch(MulVarFunc):
+    def _eval_derivative(self, s):
+        return switch(*[arg.diff(s) for arg in self.args[0:len(self.args) - 1]], self.args[-1])
+
+    def _numpycode(self, printer, **kwargs):
+        return r'SolCF.switch(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
+
+
 class Saturation(MulVarFunc):
     r"""
     The element-wise saturation a number
@@ -391,9 +382,14 @@ class Saturation(MulVarFunc):
             \end{cases}
     """
 
-    @classmethod
-    def eval(cls, v, vmin, vmax):
-        return v * In(v, vmin, vmax) + vmax * GreaterThan(v, vmax) + vmin * LessThan(v, vmin)
+    def fdiff(self, argindex=1):
+        if argindex == 1:
+            return In(*self.args)
+        else:
+            return Integer(0)
+
+    def _numpycode(self, printer, **kwargs):
+        return r'SolCF.Saturation(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 class AntiWindUp(MulVarFunc):
@@ -426,6 +422,7 @@ class AntiWindUp(MulVarFunc):
             \operatorname{AntiWindUp}(u_{\text {des }}(t), u_{\min }, u_{\max }, e(t)).
 
     """
+    arglength = 4
 
     @classmethod
     def eval(cls, u, umin, umax, e):
@@ -449,6 +446,7 @@ class Min(MulVarFunc):
         \end{cases}\end{split}
 
     """
+    arglength = 2
 
     @classmethod
     def eval(cls, x, y):
@@ -460,6 +458,7 @@ class In(MulVarFunc):
     In(v, vmin, vmax)
         return True if vmin<=v<=vmax
     """
+    arglength = 3
 
     def _eval_derivative(self, s):
         return Integer(0)
@@ -470,19 +469,14 @@ def _sympystr(self, printer, **kwargs):
                                                     op3=printer._print(self.args[2]))
 
     def _numpycode(self, printer, **kwargs):
-        return r'SolIn(' + ', '.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)
+        return r'SolCF.In(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 class GreaterThan(MulVarFunc):
     """
     Represents > operator
     """
+    arglength = 2
 
     def _eval_derivative(self, s):
         return Integer(0)
@@ -492,19 +486,14 @@ def _sympystr(self, printer, **kwargs):
                                           op2=printer._print(self.args[1]))
 
     def _numpycode(self, printer, **kwargs):
-        return r'SolGreaterThan(' + ', '.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)
+        return r'SolCF.GreaterThan(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 class LessThan(MulVarFunc):
     """
     Represents < operator
     """
+    arglength = 2
 
     def _eval_derivative(self, s):
         return Integer(0)
@@ -514,19 +503,14 @@ def _sympystr(self, printer, **kwargs):
                                           op2=printer._print(self.args[1]))
 
     def _numpycode(self, printer, **kwargs):
-        return r'SolLessThan(' + ', '.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)
+        return r'SolCF.LessThan(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 class And(MulVarFunc):
     """
     Represents bitwise_and
     """
+    arglength = 2
 
     def _eval_derivative(self, s):
         return Integer(0)
@@ -536,19 +520,14 @@ def _sympystr(self, printer, **kwargs):
                                           op2=printer._print(self.args[1]))
 
     def _numpycode(self, printer, **kwargs):
-        return r'And(' + ', '.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)
+        return r'SolCF.And(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 class Or(MulVarFunc):
     """
     Represents bitwise_or
     """
+    arglength = 2
 
     def _eval_derivative(self, s):
         return Integer(0)
@@ -558,34 +537,7 @@ def _sympystr(self, printer, **kwargs):
                                           op2=printer._print(self.args[1]))
 
     def _numpycode(self, printer, **kwargs):
-        return r'Or(' + ', '.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)
-
-
-class Not(MulVarFunc):
-    """
-    Represents bitwise_not
-    """
-
-    def _eval_derivative(self, s):
-        return Integer(0)
-
-    def _sympystr(self, printer, **kwargs):
-        return '({op1})'.format(op1=printer._print(self.args[0]))
-
-    def _numpycode(self, printer, **kwargs):
-        return r'Not(' + ', '.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)
+        return r'SolCF.Or(' + ', '.join([printer._print(arg, **kwargs) for arg in self.args]) + r')'
 
 
 # %% custom func of equation printer
@@ -605,7 +557,7 @@ def eval(cls, *args):
             raise TypeError(f"CSC_array takes 1 positional arguments but {len(args)} were given!")
 
     def _numpycode(self, printer, **kwargs):
-        return r'csc_array(' + printer._print(self.args[0]) + r')'
+        return r'scipy.sparse.csc_array(' + printer._print(self.args[0]) + r')'
 
     def _pythoncode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
@@ -624,7 +576,73 @@ def eval(cls, *args):
             raise ValueError(f"Solverz' arange object takes 2 positional arguments but {len(args)} were given!")
 
     def _numpycode(self, printer, **kwargs):
-        return r'arange(' + ', '.join([printer._print(arg) for arg in self.args]) + r')'
+        return r'np.arange(' + ', '.join([printer._print(arg) for arg in self.args]) + r')'
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+
+class Ones(Function):
+    r"""
+    The all-one vector to broadcast scalars to a vector
+    For example, 2 -> 2*Ones(eqn_size)
+    The derivative is d(Ones(eqn_size))/dx=0
+    """
+
+    @classmethod
+    def eval(cls, x):
+        if not isinstance(x, Integer):
+            raise ValueError("The arg of Ones() should be integer.")
+
+    def _eval_derivative(self, s):
+        return 0
+
+    def _numpycode(self, printer, **kwargs):
+        x = self.args[0]
+        return r'np.ones(' + printer._print(x, **kwargs) + r')'
+
+    def _lambdacode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+
+class coo_array(Function):
+
+    @classmethod
+    def eval(cls, *args):
+        if len(args) > 1:
+            raise ValueError(
+                f"Solverz' coo_array object accepts only one inputs.")
+
+    def _numpycode(self, printer, **kwargs):
+        return f'sps.coo_array({printer._print(self.args[0])})'
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+
+class extend(Function):
+
+    def _numpycode(self, printer, **kwargs):
+        return f'{printer._print(self.args[0])}.extend({printer._print(self.args[1])})'
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+
+class zeros(Function):
+    # print zeros(6,6) as zeros((6,6))
+    # or zeros(6,) as zeros((6,))
+    def _numpycode(self, printer, **kwargs):
+        if len(self.args) == 2:
+            temp1 = printer._print(self.args[0])
+            temp2 = printer._print(self.args[1])
+            return r'np.zeros((' + temp1 + ', ' + temp2 + r'))'
+        elif len(self.args) == 1:
+            temp = printer._print(self.args[0])
+            return r'np.zeros((' + temp + ', ' + r'))'
 
     def _pythoncode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
diff --git a/Solverz/sym_algebra/matrix_calculus.py b/Solverz/sym_algebra/matrix_calculus.py
index 1329a06..4d4fde9 100644
--- a/Solverz/sym_algebra/matrix_calculus.py
+++ b/Solverz/sym_algebra/matrix_calculus.py
@@ -9,6 +9,7 @@
 
 from Solverz.sym_algebra.symbols import IdxVar, IdxPara, iVar, Para
 from Solverz.sym_algebra.functions import Mat_Mul, Diag, transpose, Abs, Sign
+from Solverz.num_api.module_parser import modules
 
 
 class TMatrix:
@@ -245,7 +246,7 @@ def obtain_dim(expr) -> int:
                 symbol_dict[symbol] = np.ones((2, 2))
             elif symbol.dim == 1:
                 symbol_dict[symbol] = np.ones((2, 1))
-    temp_expr = sp.lambdify(symbol_dict.keys(), expr)
+    temp_expr = sp.lambdify(symbol_dict.keys(), expr, modules=modules)
     return temp_expr(*symbol_dict.values()).shape[1]
 
 
diff --git a/Solverz/sym_algebra/symbols.py b/Solverz/sym_algebra/symbols.py
index 37bf5cd..ed3ae8a 100644
--- a/Solverz/sym_algebra/symbols.py
+++ b/Solverz/sym_algebra/symbols.py
@@ -246,3 +246,33 @@ def _lambdacode(self, printer, **kwargs):
 
     def _pythoncode(self, printer, **kwargs):
         return self._numpycode(printer, **kwargs)
+
+
+class coo_2_csc(Symbol):
+
+    def __new__(cls, eqn_size: int, vsize: int):
+        obj = Symbol.__new__(cls, f'coo_2_csc')
+        obj.eqn_size = eqn_size
+        obj.vsize = vsize
+        return obj
+
+    def _numpycode(self, printer, **kwargs):
+        return f'sps.coo_array((data, (row, col)), ({self.eqn_size}, {self.vsize})).tocsc()'
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
+
+
+class coo_2_csc_hvp(Symbol):
+
+    def __new__(cls, eqn_size: int, vsize: int):
+        obj = Symbol.__new__(cls, f'coo_2_csc')
+        obj.eqn_size = eqn_size
+        obj.vsize = vsize
+        return obj
+
+    def _numpycode(self, printer, **kwargs):
+        return f'sps.coo_array((data_hvp, (row_hvp, col_hvp)), ({self.eqn_size}, {self.vsize})).tocsc()'
+
+    def _pythoncode(self, printer, **kwargs):
+        return self._numpycode(printer, **kwargs)
diff --git a/Solverz/sym_algebra/test/test_finite_difference.py b/Solverz/sym_algebra/test/test_finite_difference.py
deleted file mode 100644
index 691d618..0000000
--- a/Solverz/sym_algebra/test/test_finite_difference.py
+++ /dev/null
@@ -1,103 +0,0 @@
-from sympy import pycode, Integer
-
-from Solverz.sym_algebra.transform import finite_difference, semi_descritize
-from Solverz.sym_algebra.symbols import iVar, Para
-
-pi = iVar('pi')
-qi = iVar('qi')
-S = Para('S')
-va = Para('va')
-D = Para('D')
-lam = Para('lam')
-
-
-# finite difference of hyperbolic PDE, the gas equation
-def test_fd():
-    fde = finite_difference(qi,
-                            S * pi,
-                            -lam * va ** 2 * qi ** 2 / (2 * D * S * pi),
-                            [pi, qi],
-                            10,
-                            'central diff')
-    assert pycode(
-        fde) == 'S*dt*(-pi_tag_0[0:10] + pi_tag_0[1:11] - pi[0:10] + pi[1:11]) + dx*(-qi_tag_0[0:10] - qi_tag_0[1:11] + qi[0:10] + qi[1:11]) + (1/4)*dt*dx*lam*va**2*(qi_tag_0[0:10] + qi_tag_0[1:11] + qi[0:10] + qi[1:11])**2/(D*S*(pi_tag_0[0:10] + pi_tag_0[1:11] + pi[0:10] + pi[1:11]))'
-
-    fde = finite_difference(pi,
-                            va ** 2 / S * qi,
-                            Integer(0),
-                            [pi, qi],
-                            10,
-                            'central diff')
-    assert pycode(
-        fde) == 'dx*(-pi_tag_0[0:10] - pi_tag_0[1:11] + pi[0:10] + pi[1:11]) + dt*va**2*(-qi_tag_0[0:10] + qi_tag_0[1:11] - qi[0:10] + qi[1:11])/S'
-
-    fde = finite_difference(qi,
-                            S * pi,
-                            -lam * va ** 2 * qi ** 2 / (2 * D * S * pi),
-                            [pi, qi],
-                            10,
-                            'euler',
-                            direction=1)
-    assert pycode(
-        fde) == 'S*dt*(-pi[0:10] + pi[1:11]) + dx*(-qi_tag_0[1:11] + qi[1:11]) + (1/2)*qi[1:11]**2*dt*dx*lam*va**2/(pi[1:11]*D*S)'
-
-    fde = finite_difference(pi,
-                            va ** 2 / S * qi,
-                            Integer(0),
-                            [pi, qi],
-                            10,
-                            'euler',
-                            direction=1)
-    assert pycode(fde) == 'dx*(-pi_tag_0[1:11] + pi[1:11]) + dt*va**2*(-qi[0:10] + qi[1:11])/S'
-
-    fde = finite_difference(qi,
-                            S * pi,
-                            -lam * va ** 2 * qi ** 2 / (2 * D * S * pi),
-                            [pi, qi],
-                            10,
-                            'euler',
-                            direction=-1)
-    assert pycode(
-        fde) == 'S*dt*(-pi[0:10] + pi[1:11]) + dx*(-qi_tag_0[0:10] + qi[0:10]) + (1/2)*qi[0:10]**2*dt*dx*lam*va**2/(pi[0:10]*D*S)'
-
-    fde = finite_difference(pi,
-                            va ** 2 / S * qi,
-                            Integer(0),
-                            [pi, qi],
-                            10,
-                            'euler',
-                            direction=-1)
-    assert pycode(fde) == 'dx*(-pi_tag_0[0:10] + pi[0:10]) + dt*va**2*(-qi[0:10] + qi[1:11])/S'
-
-
-def test_sd():
-    # semi-discretize
-    eqn_dict = semi_descritize(qi,
-                               S * pi,
-                               -lam * va ** 2 * qi ** 2 / (2 * D * S * pi),
-                               [pi, qi],
-                               11,
-                               'TVD1',
-                               a0=va,
-                               a1=va)
-
-    rhs = eqn_dict['Ode'][0]
-    diff_var = eqn_dict['Ode'][1]
-
-    assert pycode(
-        rhs) == '-1/2*S*(-pi[0:10] + pi[2:12])/dx + (1/2)*va*(qi[0:10] - 2*qi[1:11] + qi[2:12])/dx - 1/2*qi[1:11]**2*lam*va**2/(pi[1:11]*D*S)'
-    assert pycode(diff_var) == 'qi[1:11]'
-
-    eqn_dict = semi_descritize(pi,
-                               va ** 2 / S * qi,
-                               Integer(0),
-                               [pi, qi],
-                               11,
-                               'TVD1',
-                               a0=va,
-                               a1=va)
-    rhs = eqn_dict['Ode'][0]
-    diff_var = eqn_dict['Ode'][1]
-
-    assert pycode(rhs) == '(1/2)*va*(pi[0:10] - 2*pi[1:11] + pi[2:12])/dx - 1/2*va**2*(-qi[0:10] + qi[2:12])/(S*dx)'
-    assert pycode(diff_var) == 'pi[1:11]'
diff --git a/Solverz/sym_algebra/test/test_functions.py b/Solverz/sym_algebra/test/test_functions.py
new file mode 100644
index 0000000..b6d17a3
--- /dev/null
+++ b/Solverz/sym_algebra/test/test_functions.py
@@ -0,0 +1,7 @@
+from ..functions import sin, cos
+from Solverz import Var
+
+def test_function_arg_parser():
+    # The functions parse the Solverz.ssymbol as symbols.
+    assert sin(Var('x')).args[0].name == 'x'
+    assert cos(Var('x')).args[0].name == 'x'
diff --git a/Solverz/sym_algebra/transform.py b/Solverz/sym_algebra/transform.py
deleted file mode 100644
index 0466d2a..0000000
--- a/Solverz/sym_algebra/transform.py
+++ /dev/null
@@ -1,172 +0,0 @@
-from sympy import simplify
-
-from Solverz.sym_algebra.symbols import Para, iAliasVar, IdxVar, iVar
-from Solverz.sym_algebra.functions import switch
-from Solverz.utilities.type_checker import is_number
-
-
-def finite_difference(diff_var, flux, source, two_dim_var, M, scheme='central diff', direction=None, dx=None):
-    M_ = M + 1  # for pretty printer of M in slice
-    if dx is None:
-        dx = Para('dx')
-    else:
-        if is_number(dx):
-            dx = dx
-        else:
-            raise TypeError(f'Input dx is not number!')
-
-    dt = Para('dt')
-    u = diff_var
-    u0 = iAliasVar(u.name + '_tag_0')
-
-    if scheme == 'central diff':
-        fui1j1 = flux.subs([(a, a[1:M_]) for a in two_dim_var])
-        fuij1 = flux.subs([(a, a[0:M_ - 1]) for a in two_dim_var])
-        fui1j = flux.subs([(a, iAliasVar(a.name + '_tag_0')[1:M_]) for a in two_dim_var])
-        fuij = flux.subs([(a, iAliasVar(a.name + '_tag_0')[0:M_ - 1]) for a in two_dim_var])
-
-        S = source.subs([(a, (a[1:M_] + a[0:M_ - 1] + iAliasVar(a.name + '_tag_0')[1:M_] + iAliasVar(a.name + '_tag_0')[
-                                                                                          0:M_ - 1]) / 4) for a in
-                         two_dim_var])
-
-        fde = dx * (u[1:M_] - u0[1:M_] + u[0:M_ - 1] - u0[0:M_ - 1]) \
-              + simplify(dt * (fui1j1 - fuij1 + fui1j - fuij)) \
-              - simplify(2 * dx * dt * S)
-
-    elif scheme == 'euler':
-        if direction == 1:
-            fui1j1 = flux.subs([(a, a[1:M_]) for a in two_dim_var])
-            fuij1 = flux.subs([(a, a[0:M_ - 1]) for a in two_dim_var])
-            S = source.subs([(a, a[1:M_]) for a in two_dim_var])
-            fde = dx * (u[1:M_] - u0[1:M_]) + simplify(dt * (fui1j1 - fuij1)) - simplify(dx * dt * S)
-        elif direction == -1:
-            fui1j1 = flux.subs([(a, a[1:M_]) for a in two_dim_var])
-            fuij1 = flux.subs([(a, a[0:M_ - 1]) for a in two_dim_var])
-            S = source.subs([(a, a[0:M_ - 1]) for a in two_dim_var])
-            fde = dx * (u[0:M_ - 1] - u0[0:M_ - 1]) + simplify(dt * (fui1j1 - fuij1)) - simplify(dx * dt * S)
-        else:
-            raise ValueError(f"Unimplemented direction {direction}!")
-    return fde
-
-
-def semi_descritize(diff_var,
-                    flux,
-                    source,
-                    two_dim_var,
-                    M,
-                    scheme='TVD1',
-                    a0=None,
-                    a1=None,
-                    dx=None):
-    M_ = M + 1  # for pretty printer of M in slice
-
-    if a0 is None:
-        a0 = Para('ajp12')
-    if a1 is None:
-        a1 = Para('ajm12')
-
-    if dx is None:
-        dx = Para('dx')
-    else:
-        if is_number(dx):
-            dx = dx
-        else:
-            raise TypeError(f'Input dx is not number!')
-
-    u = diff_var
-    if scheme == 'TVD2':
-        # j=1
-        # f(u[2])
-        fu2 = flux.subs([(var, var[2]) for var in two_dim_var])
-        # f(u[0])=f(2*uL-u[1])
-        fu0 = flux.subs([(var, var[0]) for var in two_dim_var])
-        # S(u[1])
-        Su1 = source.subs([(var, var[1]) for var in two_dim_var])
-        ode_rhs1 = -simplify((fu2 - fu0) / (2 * dx)) \
-                   + simplify((a0[0] * (u[2] - u[1]) - a1[0] * (u[1] - u[0])) / (2 * dx)) \
-                   + simplify(Su1)
-
-        # j=M-1
-        # f(u[M])=f(2*uR-u[M-1])
-        fum = flux.subs([(var, var[M]) for var in two_dim_var])
-        # f(u[M-2])
-        fum2 = flux.subs([(var, var[M - 2]) for var in two_dim_var])
-        # S(u[M-1])
-        SuM1 = source.subs([(var, var[M - 1]) for var in two_dim_var])
-        ode_rhs3 = -simplify((fum - fum2) / (2 * dx)) \
-                   + simplify((a0[-1] * (u[M] - u[M - 1]) - a1[-1] * (u[M - 1] - u[M - 2])) / (2 * dx)) \
-                   + simplify(SuM1)
-
-        # 2<=j<=M-2
-        def ujprime(U: IdxVar, v: int):
-            # for given u_j,
-            # returns
-            # u^+_{j+1/2} case v==0,
-            # u^-_{j+1/2} case 1,
-            # u^+_{j-1/2} case 2,
-            # u^-_{j-1/2} case 3
-            if not isinstance(U.index, slice):
-                raise TypeError("Index of IdxVar must be slice object")
-            start = U.index.start
-            stop = U.index.stop
-            step = U.index.step
-            U = U.symbol0
-            Ux = iVar(U.name + 'x')
-
-            # u_j
-            Uj = U[start:stop:step]
-            # (u_x)_j
-            Uxj = Ux[start:stop:step]
-            # u_{j+1}
-            Ujp1 = U[start + 1:stop + 1:step]
-            # (u_x)_{j+1}
-            Uxjp1 = Ux[start + 1:stop + 1:step]
-            # u_{j-1}
-            Ujm1 = U[start - 1:stop - 1:step]
-            # (u_x)_{j-1}
-            Uxjm1 = Ux[start - 1:stop - 1:step]
-
-            if v == 0:
-                return Ujp1 - dx / 2 * Uxjp1
-            elif v == 1:
-                return Uj + dx / 2 * Uxj
-            elif v == 2:
-                return Uj - dx / 2 * Uxj
-            elif v == 3:
-                return Ujm1 + dx / 2 * Uxjm1
-            else:
-                raise ValueError("v=0 or 1 or 2 or 3!")
-
-        # j\in [2:M_-2]
-        Suj = source.subs([(var, var[2:M_ - 2]) for var in two_dim_var])
-        Hp = (flux.subs([(var, ujprime(var[2:M_ - 2], 0)) for var in two_dim_var]) +
-              flux.subs([(var, ujprime(var[2:M_ - 2], 1)) for var in two_dim_var])) / 2 \
-             - a0[2:M_ - 2] / 2 * (ujprime(u[2:M_ - 2], 0) - ujprime(u[2:M_ - 2], 1))
-        Hm = (flux.subs([(var, ujprime(var[2:M_ - 2], 2)) for var in two_dim_var]) +
-              flux.subs([(var, ujprime(var[2:M_ - 2], 3)) for var in two_dim_var])) / 2 \
-             - a1[2:M_ - 2] / 2 * (ujprime(u[2:M_ - 2], 2) - ujprime(u[2:M_ - 2], 3))
-        ode_rhs2 = -simplify(Hp - Hm) / dx + Suj
-
-        theta = Para('theta')
-        ux = iVar(u.name + 'x')
-        minmod_flag = Para('minmod_flag_of_' + ux.name)
-        minmod_rhs = ux[1:M_ - 1] - switch(theta * (u[1:M_ - 1] - u[0:M_ - 2]) / dx,
-                                           (u[2:M_] - u[0:M_ - 2]) / (2 * dx),
-                                           theta * (u[2:M_] - u[1:M_ - 1]) / dx,
-                                           0,
-                                           minmod_flag)
-
-        return {'Ode': [(ode_rhs1, u[1]), (ode_rhs2, u[2:M_ - 2]), (ode_rhs3, u[M - 1])],
-                'Eqn': [minmod_rhs, ux[0], ux[M]]}
-    elif scheme == 'TVD1':
-        # 1<=j<=M-1
-        # f(u[j+1])
-        fu1 = flux.subs([(var, var[2:M_]) for var in two_dim_var])
-        # f(u[j-1])
-        fu2 = flux.subs([(var, var[0:M_ - 2]) for var in two_dim_var])
-        # S(u[j])
-        Su = source.subs([(var, var[1:M_ - 1]) for var in two_dim_var])
-        ode_rhs = -simplify((fu1 - fu2) / (2 * dx)) \
-                  + simplify((a0 * (u[2:M_] - u[1:M_ - 1]) - a1 * (u[1:M_ - 1] - u[0:M_ - 2])) / (2 * dx)) \
-                  + simplify(Su)
-        return {'Ode': (ode_rhs, u[1:M_ - 1])}
diff --git a/Solverz/variable/ssymbol.py b/Solverz/variable/ssymbol.py
index edba549..7c7c13e 100644
--- a/Solverz/variable/ssymbol.py
+++ b/Solverz/variable/ssymbol.py
@@ -110,7 +110,7 @@ def __init__(self, name: str, value=None, init=None):
 
 class AliasVar(sSymBasic):
     def __init__(self, name: str, step=1, value=None, init=None):
-        name = name + '_tag_' + str(step-1)
+        name = name + '_tag_' + str(step - 1)
         super().__init__(name=name, Type='iAliasVar', value=value, dim=1, init=init)
         self.step = step
 
@@ -125,3 +125,9 @@ def sSym2Sym(x):
     else:
         raise TypeError(f'Input type {type(x)} not supported!')
 
+
+def sVar2Var(var: Union[Var, iVar, List[Union[iVar, Var]]]) -> Union[iVar, List[iVar]]:
+    if isinstance(var, list):
+        return [arg.symbol if isinstance(arg, Var) else arg for arg in var]
+    else:
+        return var.symbol if isinstance(var, Var) else var
diff --git a/docs/src/gettingstart.md b/docs/src/gettingstart.md
index dc3d6bb..a2acf10 100644
--- a/docs/src/gettingstart.md
+++ b/docs/src/gettingstart.md
@@ -6,7 +6,7 @@
 Solverz supports three types of equality constraints, which are respectively
 
 1. Algebraic Equations (AEs) with general formula $0=F(y,p)$
-2. Finite Difference Algebraic Equations (FDAEs) with general formula $0=F(y,p,y_0)$
+2. Finite Difference Algebraic Equations (FDAEs) with general formula $0=F(t,y,p,y_0)$
 3. Differential Algebraic Equations (DAEs) with general formula $M\dot{y}=F(t,y,p)$
 
 where $y$ is the known variable, $F$ is the mapping, $p$ is the parameter set of your models, $y_0$ is the previous time node value of $y$, $M$is the singular mass matrix.
diff --git a/docs/src/intro.md b/docs/src/intro.md
index dafb9e1..83ae58c 100644
--- a/docs/src/intro.md
+++ b/docs/src/intro.md
@@ -4,14 +4,6 @@
 
 Solverz aims to help you model and solve your equations more efficiently.
 
-Solverz supports three types of abstract equation types, that are
-
-- Algebraic Equations (AEs) $0=F(y,p)$
-- Finite Difference Algebraic Equations (FDAEs) $0=F(y,p,y_0)$
-- Differential Algebraic Equations (DAEs) $M\dot{y}=F(t,y,p)$
-
-where $p$ is the parameter set of your models, $y_0$ is the previous time node value of $y$.
-
 Say, we want to know how long it takes for an apple to fall from a tree to the ground. We have the 
 differential equations 
 
diff --git a/pyproject.toml b/pyproject.toml
index 4713733..be75190 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -36,6 +36,9 @@ classifiers = [
     "Operating System :: OS Independent",
 ]
 
+[tool.pytest.ini_options]
+addopts = "-vv"
+
 [project.urls]
 Homepage = "https://github.com/smallbunnies/Solverz"
 Issues = "https://github.com/smallbunnies/Solverz/issues"