diff --git a/Solverz/__init__.py b/Solverz/__init__.py
index f737da7..ef7cf9f 100644
--- a/Solverz/__init__.py
+++ b/Solverz/__init__.py
@@ -3,11 +3,11 @@
 from Solverz.equation.param import Param, IdxParam, TimeSeriesParam
 from Solverz.sym_algebra.symbols import idx, Para, Var, AliasVar
 from Solverz.sym_algebra.functions import Sign, Abs, transpose, exp, Diag, Mat_Mul, sin, cos
-from Solverz.numerical_interface.custom_function import minmod_flag, minmod
+from Solverz.num_api.custom_function import minmod_flag, minmod
 from Solverz.variable.variables import Vars, TimeVars, as_Vars
 from Solverz.solvers.nlaesolver import nr_method, continuous_nr
 from Solverz.solvers.daesolver import Rodas, Opt, implicit_trapezoid, backward_euler
 from Solverz.solvers.fdesolver import fdae_solver, fdae_ss_solver
-from Solverz.numerical_interface.num_eqn import made_numerical, parse_dae_v, parse_ae_v, render_modules
+from Solverz.code_printer.make_pyfunc import made_numerical
 from Solverz.utilities.io import save, load
 from Solverz.utilities.profile import count_time
diff --git a/Solverz/numerical_interface/__init__.py b/Solverz/code_printer/__init__.py
similarity index 100%
rename from Solverz/numerical_interface/__init__.py
rename to Solverz/code_printer/__init__.py
diff --git a/Solverz/code_printer/make_module.py b/Solverz/code_printer/make_module.py
new file mode 100644
index 0000000..a9c685b
--- /dev/null
+++ b/Solverz/code_printer/make_module.py
@@ -0,0 +1,32 @@
+from typing import List
+
+from Solverz.code_printer.py_printer import render_modules as py_render
+from Solverz.equation.equations import AE, FDAE, DAE
+from Solverz.variable.variables import Vars
+
+
+class module_printer:
+    def __init__(self,
+                 mdl: AE | FDAE | DAE,
+                 variables: Vars | List[Vars],
+                 name: str,
+                 lang='python',
+                 directory=None,
+                 jit=False):
+        self.name = name
+        self.lang = lang
+        self.mdl = mdl
+        if isinstance(variables, Vars):
+            self.variables = [variables]
+        else:
+            self.variables = variables
+        self.directory = directory
+        self.jit = jit
+
+    def render(self):
+        if self.lang == 'python':
+            py_render(self.mdl,
+                      *self.variables,
+                      name=self.name,
+                      directory=self.directory,
+                      numba=self.jit)
diff --git a/Solverz/code_printer/make_pyfunc.py b/Solverz/code_printer/make_pyfunc.py
new file mode 100644
index 0000000..4c664e7
--- /dev/null
+++ b/Solverz/code_printer/make_pyfunc.py
@@ -0,0 +1,5 @@
+import numpy as np
+
+from Solverz.code_printer.py_printer import made_numerical
+
+
diff --git a/Solverz/numerical_interface/test/__init__.py b/Solverz/code_printer/matlab_printer.py
similarity index 100%
rename from Solverz/numerical_interface/test/__init__.py
rename to Solverz/code_printer/matlab_printer.py
diff --git a/Solverz/numerical_interface/code_printer.py b/Solverz/code_printer/py_printer.py
similarity index 87%
rename from Solverz/numerical_interface/code_printer.py
rename to Solverz/code_printer/py_printer.py
index bfa5255..42ecf26 100644
--- a/Solverz/numerical_interface/code_printer.py
+++ b/Solverz/code_printer/py_printer.py
@@ -14,15 +14,125 @@
 from numbers import Number
 
 from Solverz.equation.equations import Equations as SymEquations, FDAE as SymFDAE, DAE as SymDAE, AE as SymAE
-from Solverz.variable.variables import Vars
+from Solverz.variable.variables import Vars, TimeVars
 from Solverz.equation.param import TimeSeriesParam
 from Solverz.sym_algebra.symbols import Var, SolDict, Para, idx, IdxSymBasic
 from Solverz.sym_algebra.functions import zeros, CSC_array, Arange
 from Solverz.utilities.address import Address
 from Solverz.utilities.io import save
-
+from Solverz.variable.variables import combine_Vars
+from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.num_eqn import nAE, nFDAE, nDAE
 
 # %%
+
+
+def parse_p(ae: SymEquations):
+    p = dict()
+    for param_name, param in ae.PARAM.items():
+        if isinstance(param, TimeSeriesParam):
+            p.update({param_name: param})
+        else:
+            p.update({param_name: param.v})
+    return p
+
+
+def parse_trigger_fun(ae: SymEquations):
+    func = dict()
+    for para_name, param in ae.PARAM.items():
+        func.update({para_name + '_trigger_func': param.trigger_fun})
+
+    return func
+
+
+def made_numerical(eqn: SymEquations, *xys, sparse=False, output_code=False):
+    """
+    factory method of numerical equations
+    """
+    print(f"Printing numerical codes of {eqn.name}")
+    eqn.assign_eqn_var_address(*xys)
+    code_F = print_F(eqn)
+    code_J = print_J(eqn, sparse)
+    custom_func = dict()
+    custom_func.update(numerical_interface)
+    custom_func.update(parse_trigger_fun(eqn))
+    F = Solverzlambdify(code_F, 'F_', modules=[custom_func, 'numpy'])
+    J = Solverzlambdify(code_J, 'J_', modules=[custom_func, 'numpy'])
+    p = parse_p(eqn)
+    print('Complete!')
+    if isinstance(eqn, SymAE) and not isinstance(eqn, SymFDAE):
+        num_eqn = nAE(F, J, p)
+    elif isinstance(eqn, SymFDAE):
+        num_eqn = nFDAE(F, J, p, eqn.nstep)
+    elif isinstance(eqn, SymDAE):
+        num_eqn = nDAE(eqn.M, F, J, p)
+    else:
+        raise ValueError(f'Unknown equation type {type(eqn)}')
+    if output_code:
+        return num_eqn, {'F': code_F, 'J': code_J}
+    else:
+        return num_eqn
+
+
+def render_modules(eqn: SymEquations, *xys, name, directory=None, numba=False):
+    """
+    factory method of numerical equations
+    """
+    print(f"Printing python codes of {eqn.name}...")
+    eqn.assign_eqn_var_address(*xys)
+    p = parse_p(eqn)
+    code_F = print_F_numba(eqn)
+    code_inner_F = print_inner_F(eqn)
+    code_sub_inner_F = print_sub_inner_F(eqn)
+    code_J = print_J_numba(eqn)
+    codes = print_inner_J(eqn, *xys)
+    code_inner_J = codes['code_inner_J']
+    code_sub_inner_J = codes['code_sub_inner_J']
+    custom_func = dict()
+    custom_func.update(numerical_interface)
+    custom_func.update(parse_trigger_fun(eqn))
+
+    print('Complete!')
+
+    eqn_parameter = {}
+    if isinstance(eqn, SymAE) and not isinstance(eqn, SymFDAE):
+        eqn_type = 'AE'
+    elif isinstance(eqn, SymFDAE):
+        eqn_type = 'FDAE'
+        eqn_parameter.update({'nstep': eqn.nstep})
+    elif isinstance(eqn, SymDAE):
+        eqn_type = 'DAE'
+        eqn_parameter.update({'M': eqn.M})
+    else:
+        raise ValueError(f'Unknown equation type {type(eqn)}')
+
+    if len(xys) == 1:
+        y = xys[0]
+    else:
+        y = xys[0]
+        for arg in xys[1:]:
+            y = combine_Vars(y, arg)
+
+    code_dict = {'F': code_F,
+                 'inner_F': code_inner_F,
+                 'sub_inner_F': code_sub_inner_F,
+                 'J': code_J,
+                 'inner_J': code_inner_J,
+                 'sub_inner_J': code_sub_inner_J}
+    eqn_parameter.update({'row': codes['row'], 'col': codes['col'], 'data': codes['data']})
+    print(f"Rendering python modules!")
+    render_as_modules(name,
+                      code_dict,
+                      eqn_type,
+                      p,
+                      eqn_parameter,
+                      y,
+                      [custom_func, 'numpy'],
+                      numba,
+                      directory)
+    print('Complete!')
+
+
 def is_valid_python_module_name(module_name):
     pattern = re.compile(r'^[a-zA-Z_][a-zA-Z0-9_]*$')
     return bool(pattern.match(module_name))
@@ -58,7 +168,7 @@ def print_init_code(eqn_type: str, module_name, eqn_param):
     code += 'import time\n'
     match eqn_type:
         case 'AE':
-            code += 'from Solverz.numerical_interface.num_eqn import nAE\n'
+            code += 'from Solverz.num_api.num_eqn import nAE\n'
             code += 'mdl = nAE(F_, J_, p)\n'
             code_compile = 'mdl.F(y.array, p)\nmdl.J(y.array, p)\n'
         case 'FDAE':
@@ -66,12 +176,12 @@ def print_init_code(eqn_type: str, module_name, eqn_param):
                 nstep = eqn_param['nstep']
             except KeyError as e:
                 raise ValueError("Cannot parse nstep attribute for FDAE object printing!")
-            code += 'from Solverz.numerical_interface.num_eqn import nFDAE\n'
+            code += 'from Solverz.num_api.num_eqn import nFDAE\n'
             code += 'mdl = nFDAE(F_, J_, p, setting["nstep"])\n'
             args_str = ', '.join(['y.array' for i in range(nstep)])
             code_compile = f'mdl.F(0, y.array, p, {args_str})\nmdl.J(0, y.array, p, {args_str})\n'
         case 'DAE':
-            code += 'from Solverz.numerical_interface.num_eqn import nDAE\n'
+            code += 'from Solverz.num_api.num_eqn import nDAE\n'
             code += 'mdl = nDAE(setting["M"], F_, J_, p)\n'
             code_compile = 'mdl.F(0, y.array, p)\nmdl.J(0, y.array, p)\n'
         case _:
@@ -516,7 +626,8 @@ def print_J_numba(ae: SymEquations):
     body.extend(param_assignments)
     body.extend(print_trigger(ae))
     body.extend([Assignment(Var('data', internal_use=True),
-                            FunctionCall('inner_J', [symbols('_data_', real=True)] + [arg.name for arg in var_list + param_list]))])
+                            FunctionCall('inner_J', [symbols('_data_', real=True)] + [arg.name for arg in
+                                                                                      var_list + param_list]))])
     body.extend([Return(coo_2_csc(ae))])
     fd = FunctionDefinition.from_FunctionPrototype(fp, body)
     return pycode(fd, fully_qualified_modules=False)
@@ -577,7 +688,8 @@ def print_F_numba(ae: SymEquations):
     param_assignments, param_list = print_param(ae, numba_printer=True)
     body.extend(param_assignments)
     body.extend(print_trigger(ae))
-    body.extend([Return(FunctionCall('inner_F', [symbols('_F_', real=True)] + [arg.name for arg in var_list + param_list]))])
+    body.extend(
+        [Return(FunctionCall('inner_F', [symbols('_F_', real=True)] + [arg.name for arg in var_list + param_list]))])
     fd = FunctionDefinition.from_FunctionPrototype(fp, body)
     return pycode(fd, fully_qualified_modules=False)
 
diff --git a/Solverz/solvers/test/test_numjac.py b/Solverz/code_printer/test/__init__.py
similarity index 100%
rename from Solverz/solvers/test/test_numjac.py
rename to Solverz/code_printer/test/__init__.py
diff --git a/Solverz/equation/eqn.py b/Solverz/equation/eqn.py
index bf73632..7cde6bb 100644
--- a/Solverz/equation/eqn.py
+++ b/Solverz/equation/eqn.py
@@ -11,7 +11,7 @@
 from Solverz.sym_algebra.functions import Mat_Mul, Slice, F
 from Solverz.sym_algebra.matrix_calculus import MixedEquationDiff
 from Solverz.sym_algebra.transform import finite_difference, semi_descritize
-from Solverz.numerical_interface.custom_function import numerical_interface
+from Solverz.num_api.custom_function import numerical_interface
 
 
 class Eqn:
diff --git a/Solverz/equation/param.py b/Solverz/equation/param.py
index 08056cf..3394f0c 100644
--- a/Solverz/equation/param.py
+++ b/Solverz/equation/param.py
@@ -4,7 +4,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-from Solverz.numerical_interface.Array import Array
+from Solverz.num_api.Array import Array
 
 
 class Param:
diff --git a/Solverz/numerical_interface/Array.py b/Solverz/num_api/Array.py
similarity index 100%
rename from Solverz/numerical_interface/Array.py
rename to Solverz/num_api/Array.py
diff --git a/Solverz/num_api/__init__.py b/Solverz/num_api/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/Solverz/numerical_interface/custom_function.py b/Solverz/num_api/custom_function.py
similarity index 100%
rename from Solverz/numerical_interface/custom_function.py
rename to Solverz/num_api/custom_function.py
diff --git a/Solverz/num_api/num_eqn.py b/Solverz/num_api/num_eqn.py
new file mode 100644
index 0000000..5ab8134
--- /dev/null
+++ b/Solverz/num_api/num_eqn.py
@@ -0,0 +1,43 @@
+from typing import Callable, Dict
+
+import numpy as np
+
+from Solverz.equation.param import TimeSeriesParam
+from Solverz.variable.variables import Vars, TimeVars, combine_Vars
+
+
+class nAE:
+
+    def __init__(self,
+                 F: Callable,
+                 J: Callable,
+                 p: Dict):
+        self.F = F
+        self.J = J
+        self.p = p
+
+
+class nFDAE:
+
+    def __init__(self,
+                 F: callable,
+                 J: callable,
+                 p: dict,
+                 nstep: int = 0):
+        self.F = F
+        self.J = J
+        self.p = p
+        self.nstep = nstep
+
+
+class nDAE:
+
+    def __init__(self,
+                 M,
+                 F: Callable,
+                 J: Callable,
+                 p: Dict):
+        self.M = M
+        self.F = F
+        self.J = J
+        self.p = p
diff --git a/Solverz/solvers/numjac.py b/Solverz/num_api/numjac.py
similarity index 61%
rename from Solverz/solvers/numjac.py
rename to Solverz/num_api/numjac.py
index 3ada6ab..32d6891 100644
--- a/Solverz/solvers/numjac.py
+++ b/Solverz/num_api/numjac.py
@@ -52,3 +52,42 @@ def numjac(F, t, y, Fty, thresh, S, g, vecon, central_diff):
     dFdyt = csc_array(df)
 
     return Fty, dFdyt, nfevals
+
+
+def numjac_ae(F, y, thresh):
+    # Add t to the end of y and adjust thresh accordingly
+
+    facmax = 0.1
+    ny = len(y)
+    fac = np.sqrt(np.finfo(float).eps) + np.zeros(ny)
+    yscale = np.maximum(0.1 * np.abs(y), thresh)
+    del_ = (y + fac * yscale) - y
+    jj = np.where(del_ == 0)[0]
+
+    for j in jj:
+        while True:
+            if fac[j] < facmax:
+                fac[j] = min(100 * fac[j], facmax)
+                del_[j] = (y[j] + fac[j] * yscale[j]) - y[j]
+                if del_[j] != 0:
+                    break
+            else:
+                del_[j] = thresh[j]
+                break
+
+    nfevals = ny
+    df = np.zeros((ny, ny))
+
+    for jj in range(0, ny):
+        ydel = y.copy()
+        ydel[jj] += del_[jj]
+
+        ydel_1 = y.copy()
+        ydel_1[jj] -= del_[jj]
+        df[:, jj] = (F(ydel) - F(ydel_1)) / (2 * del_[jj])
+        nfevals += 1
+
+    # Convert df to sparse matrix dFdyt
+    dFdyt = csc_array(df)
+
+    return dFdyt, nfevals
diff --git a/Solverz/num_api/test/__init__.py b/Solverz/num_api/test/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/Solverz/numerical_interface/test/test_Array.py b/Solverz/num_api/test/test_Array.py
similarity index 98%
rename from Solverz/numerical_interface/test/test_Array.py
rename to Solverz/num_api/test/test_Array.py
index 05bc952..6128dfd 100644
--- a/Solverz/numerical_interface/test/test_Array.py
+++ b/Solverz/num_api/test/test_Array.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.sparse import csc_array
 
-from Solverz.numerical_interface.Array import Array
+from Solverz.num_api.Array import Array
 
 
 def test_Array():
diff --git a/Solverz/numerical_interface/test/test_code_printer.py b/Solverz/num_api/test/test_code_printer.py
similarity index 98%
rename from Solverz/numerical_interface/test/test_code_printer.py
rename to Solverz/num_api/test/test_code_printer.py
index 0a6651c..7518001 100644
--- a/Solverz/numerical_interface/test/test_code_printer.py
+++ b/Solverz/num_api/test/test_code_printer.py
@@ -1,6 +1,6 @@
 from sympy import symbols, pycode, Integer
 
-from Solverz.numerical_interface.code_printer import _parse_jac_eqn_address, _parse_jac_var_address, _parse_jac_data, \
+from Solverz.code_printer.py_printer import _parse_jac_eqn_address, _parse_jac_var_address, _parse_jac_data, \
     print_J_block, _print_F_assignment, _print_var_parser
 from Solverz.sym_algebra.symbols import idx, Var, Para
 
diff --git a/Solverz/num_api/test/test_numjac.py b/Solverz/num_api/test/test_numjac.py
new file mode 100644
index 0000000..25f2858
--- /dev/null
+++ b/Solverz/num_api/test/test_numjac.py
@@ -0,0 +1,25 @@
+import numpy as np
+
+from Solverz.num_api.numjac import numjac_ae
+
+
+def f(x):
+    a = 1 + np.exp(x[0]) + np.sin(x[1])
+    b = x[1] ** 2 + np.cos(x[0])
+    return np.array([a, b])
+
+
+def Janaly(x):
+    return np.array([[np.exp(x[0]), np.cos(x[1])],
+                     [-np.sin(x[0]), 2 * x[1]]])
+
+
+def Jnum(x):
+    return numjac_ae(lambda xi: f(xi),
+                     x,
+                     np.ones(len(x)) * 1.e-12)[0]
+
+
+point = np.array([1.0, 2.0])
+
+np.testing.assert_allclose(Jnum(point).toarray(), Janaly(point), rtol=1e-5)
diff --git a/Solverz/numerical_interface/num_eqn.py b/Solverz/numerical_interface/num_eqn.py
deleted file mode 100644
index 02df32c..0000000
--- a/Solverz/numerical_interface/num_eqn.py
+++ /dev/null
@@ -1,166 +0,0 @@
-from typing import Callable, Dict, Literal
-
-import numpy as np
-
-from Solverz.equation.equations import Equations as SymEquations, AE as SymAE, DAE as SymDAE, \
-    FDAE as SymFDAE
-from Solverz.equation.param import TimeSeriesParam
-from Solverz.numerical_interface.code_printer import print_J, print_F, Solverzlambdify, print_F_numba, print_J_numba, \
-    print_inner_F, print_inner_J, render_as_modules, print_sub_inner_F
-from Solverz.numerical_interface.custom_function import numerical_interface
-from Solverz.utilities.address import Address
-from Solverz.variable.variables import Vars, TimeVars, combine_Vars
-
-
-def parse_ae_v(y: np.ndarray, var_address: Address):
-    return Vars(var_address, y)
-
-
-def parse_dae_v(y: np.ndarray, var_address: Address):
-    temp = Vars(var_address, y[0, :])
-    temp = TimeVars(temp, y.shape[0])
-    temp.array[:, :] = y
-    return temp
-
-
-class nAE:
-
-    def __init__(self,
-                 F: Callable,
-                 J: Callable,
-                 p: Dict):
-        self.F = F
-        self.J = J
-        self.p = p
-
-
-class nFDAE:
-
-    def __init__(self,
-                 F: callable,
-                 J: callable,
-                 p: dict,
-                 nstep: int = 0):
-        self.F = F
-        self.J = J
-        self.p = p
-        self.nstep = nstep
-
-
-class nDAE:
-
-    def __init__(self,
-                 M,
-                 F: Callable,
-                 J: Callable,
-                 p: Dict):
-        self.M = M
-        self.F = F
-        self.J = J
-        self.p = p
-
-
-def parse_p(ae: SymEquations):
-    p = dict()
-    for param_name, param in ae.PARAM.items():
-        if isinstance(param, TimeSeriesParam):
-            p.update({param_name: param})
-        else:
-            p.update({param_name: param.v})
-    return p
-
-
-def parse_trigger_fun(ae: SymEquations):
-    func = dict()
-    for para_name, param in ae.PARAM.items():
-        func.update({para_name + '_trigger_func': param.trigger_fun})
-
-    return func
-
-
-def made_numerical(eqn: SymEquations, *xys, sparse=False, output_code=False):
-    """
-    factory method of numerical equations
-    """
-    print(f"Printing numerical codes of {eqn.name}")
-    eqn.assign_eqn_var_address(*xys)
-    code_F = print_F(eqn)
-    code_J = print_J(eqn, sparse)
-    custom_func = dict()
-    custom_func.update(numerical_interface)
-    custom_func.update(parse_trigger_fun(eqn))
-    F = Solverzlambdify(code_F, 'F_', modules=[custom_func, 'numpy'])
-    J = Solverzlambdify(code_J, 'J_', modules=[custom_func, 'numpy'])
-    p = parse_p(eqn)
-    print('Complete!')
-    if isinstance(eqn, SymAE) and not isinstance(eqn, SymFDAE):
-        num_eqn = nAE(F, J, p)
-    elif isinstance(eqn, SymFDAE):
-        num_eqn = nFDAE(F, J, p, eqn.nstep)
-    elif isinstance(eqn, SymDAE):
-        num_eqn = nDAE(eqn.M, F, J, p)
-    else:
-        raise ValueError(f'Unknown equation type {type(eqn)}')
-    if output_code:
-        return num_eqn, {'F': code_F, 'J': code_J}
-    else:
-        return num_eqn
-
-
-def render_modules(eqn: SymEquations, *xys, name, directory=None, numba=False):
-    """
-    factory method of numerical equations
-    """
-    print(f"Printing python codes of {eqn.name}...")
-    eqn.assign_eqn_var_address(*xys)
-    p = parse_p(eqn)
-    code_F = print_F_numba(eqn)
-    code_inner_F = print_inner_F(eqn)
-    code_sub_inner_F = print_sub_inner_F(eqn)
-    code_J = print_J_numba(eqn)
-    codes = print_inner_J(eqn, *xys)
-    code_inner_J = codes['code_inner_J']
-    code_sub_inner_J = codes['code_sub_inner_J']
-    custom_func = dict()
-    custom_func.update(numerical_interface)
-    custom_func.update(parse_trigger_fun(eqn))
-
-    print('Complete!')
-
-    eqn_parameter = {}
-    if isinstance(eqn, SymAE) and not isinstance(eqn, SymFDAE):
-        eqn_type = 'AE'
-    elif isinstance(eqn, SymFDAE):
-        eqn_type = 'FDAE'
-        eqn_parameter.update({'nstep': eqn.nstep})
-    elif isinstance(eqn, SymDAE):
-        eqn_type = 'DAE'
-        eqn_parameter.update({'M': eqn.M})
-    else:
-        raise ValueError(f'Unknown equation type {type(eqn)}')
-
-    if len(xys) == 1:
-        y = xys[0]
-    else:
-        y = xys[0]
-        for arg in xys[1:]:
-            y = combine_Vars(y, arg)
-
-    code_dict = {'F': code_F,
-                 'inner_F': code_inner_F,
-                 'sub_inner_F': code_sub_inner_F,
-                 'J': code_J,
-                 'inner_J': code_inner_J,
-                 'sub_inner_J': code_sub_inner_J}
-    eqn_parameter.update({'row': codes['row'], 'col': codes['col'], 'data': codes['data']})
-    print(f"Rendering python modules!")
-    render_as_modules(name,
-                      code_dict,
-                      eqn_type,
-                      p,
-                      eqn_parameter,
-                      y,
-                      [custom_func, 'numpy'],
-                      numba,
-                      directory)
-    print('Complete!')
diff --git a/Solverz/solvers/daesolver.py b/Solverz/solvers/daesolver.py
index 53c7a3c..a4c15e3 100644
--- a/Solverz/solvers/daesolver.py
+++ b/Solverz/solvers/daesolver.py
@@ -11,7 +11,7 @@
 from Solverz.solvers.nlaesolver import nr_method
 from Solverz.sym_algebra.symbols import Var
 from Solverz.variable.variables import TimeVars, Vars, as_Vars, combine_Vars
-from Solverz.numerical_interface.num_eqn import nDAE, nAE
+from Solverz.num_api.num_eqn import nDAE, nAE
 from Solverz.solvers.stats import Stats
 from Solverz.solvers.option import Opt
 from Solverz.solvers.parser import dae_io_parser
diff --git a/Solverz/solvers/fdesolver.py b/Solverz/solvers/fdesolver.py
index b2c92d3..47b21af 100644
--- a/Solverz/solvers/fdesolver.py
+++ b/Solverz/solvers/fdesolver.py
@@ -3,7 +3,7 @@
 import numpy as np
 import tqdm
 
-from Solverz.numerical_interface.num_eqn import nFDAE, nAE
+from Solverz.num_api.num_eqn import nFDAE, nAE
 from Solverz.solvers.nlaesolver import nr_method
 from Solverz.solvers.stats import Stats
 from Solverz.solvers.option import Opt
diff --git a/Solverz/solvers/nlaesolver.py b/Solverz/solvers/nlaesolver.py
index f7d74fb..5ec264d 100644
--- a/Solverz/solvers/nlaesolver.py
+++ b/Solverz/solvers/nlaesolver.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numpy import abs, max
 
-from Solverz.numerical_interface.num_eqn import nAE
+from Solverz.num_api.num_eqn import nAE
 from Solverz.solvers.laesolver import solve
 from Solverz.variable.variables import Vars
 from Solverz.solvers.option import Opt
diff --git a/Solverz/solvers/parser.py b/Solverz/solvers/parser.py
index 9784dec..8364b75 100644
--- a/Solverz/solvers/parser.py
+++ b/Solverz/solvers/parser.py
@@ -1,7 +1,8 @@
 from typing import Callable, Any, List
 import numpy as np
-from Solverz.variable.variables import Vars
-from Solverz.numerical_interface.num_eqn import nAE, nDAE, nFDAE, parse_dae_v
+from Solverz.variable.variables import Vars, TimeVars
+from Solverz.utilities.address import Address
+from Solverz.num_api.num_eqn import nAE, nDAE, nFDAE
 from Solverz.solvers.option import Opt
 
 
@@ -32,7 +33,7 @@ def wrapper(eqn: nAE, y0: np.ndarray | Vars, opt: Opt = None):
 
         # Wrap the output in Vars if the original y0 was a Vars instance
         if original_y0_is_vars:
-            y = Vars(y0.a, y)  # Assume y0.a is accessible and relevant
+            y = parse_ae_v(y, y0.a)  # Assume y0.a is accessible and relevant
 
         # Return results, with stats if they were provided
         return (y, stats) if stats is not None else y
@@ -94,3 +95,14 @@ def wrapper(eqn: nDAE, tspan: List | np.ndarray, y0: np.ndarray | Vars, opt: Opt
         return T, y, stats
 
     return wrapper
+
+
+def parse_ae_v(y: np.ndarray, var_address: Address):
+    return Vars(var_address, y)
+
+
+def parse_dae_v(y: np.ndarray, var_address: Address):
+    temp = Vars(var_address, y[0, :])
+    temp = TimeVars(temp, y.shape[0])
+    temp.array[:, :] = y
+    return temp
diff --git a/Solverz/sym_algebra/symbols.py b/Solverz/sym_algebra/symbols.py
index 29f1baf..3167d7e 100644
--- a/Solverz/sym_algebra/symbols.py
+++ b/Solverz/sym_algebra/symbols.py
@@ -3,7 +3,7 @@
 import numpy as np
 from sympy import Symbol, Expr
 
-from Solverz.numerical_interface.Array import Array
+from Solverz.num_api.Array import Array
 
 Sympify_Mapping = {}
 
diff --git a/Solverz/sym_algebra/test/test_num_alg.py b/Solverz/sym_algebra/test/test_num_alg.py
index 73452c8..ef74cd8 100644
--- a/Solverz/sym_algebra/test/test_num_alg.py
+++ b/Solverz/sym_algebra/test/test_num_alg.py
@@ -3,7 +3,7 @@
 import inspect
 
 from Solverz.sym_algebra.symbols import Var, Para, idx, IdxVar, IdxPara, Idxidx
-from Solverz.numerical_interface.custom_function import sol_slice
+from Solverz.num_api.custom_function import sol_slice
 
 
 # test of IndexPrinter
diff --git a/Solverz/variable/variables.py b/Solverz/variable/variables.py
index a5dab43..f39d636 100644
--- a/Solverz/variable/variables.py
+++ b/Solverz/variable/variables.py
@@ -8,7 +8,7 @@
 
 from Solverz.sym_algebra.symbols import Var
 from Solverz.utilities.address import Address, combine_Address
-from Solverz.numerical_interface.Array import Array
+from Solverz.num_api.Array import Array
 
 
 def as_Vars(var: Union[Var, List[Var]]) -> Vars:
diff --git a/tests/test_0.py b/tests/test_0.py
index 45e3a35..ae59bce 100644
--- a/tests/test_0.py
+++ b/tests/test_0.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from Solverz import Eqn, AE, nr_method, as_Vars, Var, made_numerical, parse_ae_v
+from Solverz import Eqn, AE, nr_method, as_Vars, Var, made_numerical
 
 x = Var(name='x', value=2)
 
diff --git a/tests/test_dhspf.py b/tests/test_dhspf.py
index 5a35009..aad2f63 100644
--- a/tests/test_dhspf.py
+++ b/tests/test_dhspf.py
@@ -1,7 +1,7 @@
 import pandas as pd
 import numpy as np
 
-from Solverz import Eqn, AE, nr_method, as_Vars, Var, Para, idx, Abs, exp, Mat_Mul, made_numerical, parse_ae_v
+from Solverz import Eqn, AE, nr_method, as_Vars, Var, Para, idx, Abs, exp, Mat_Mul, made_numerical
 
 # %% initialize variables and params
 sys_df = pd.read_excel('instances/4node3pipe.xlsx',
diff --git a/tests/test_dhspf_change_sign.py b/tests/test_dhspf_change_sign.py
index 7dccaff..fded905 100644
--- a/tests/test_dhspf_change_sign.py
+++ b/tests/test_dhspf_change_sign.py
@@ -3,7 +3,7 @@
 from functools import partial
 
 from Solverz import Eqn, AE, nr_method, as_Vars, Var, Param, continuous_nr, Para, idx, Abs, exp, Mat_Mul, \
-    made_numerical, parse_ae_v
+    made_numerical
 
 # %% initialize variables and params
 sys_df = pd.read_excel('instances/4node3pipe_change_sign.xlsx',
diff --git a/tests/test_dyn_gas_flow_single_pipe.py b/tests/test_dyn_gas_flow_single_pipe.py
index 29c2a50..8e3bc59 100644
--- a/tests/test_dyn_gas_flow_single_pipe.py
+++ b/tests/test_dyn_gas_flow_single_pipe.py
@@ -1,8 +1,7 @@
 import numpy as np
 import pandas as pd
 
-from Solverz import idx, Var, Para, as_Vars, HyperbolicPde, fdae_solver, Eqn, Param, IdxParam, FDAE, made_numerical, \
-    parse_dae_v, Opt
+from Solverz import idx, Var, Para, as_Vars, HyperbolicPde, fdae_solver, Eqn, Param, IdxParam, FDAE, made_numerical, Opt
 
 results = pd.read_excel('instances/dynamic_gas_flow_single_pipe.xlsx',
                         sheet_name=None,
@@ -44,8 +43,7 @@
 gas_FDE.param_initializer('dt', Param('dt', value=5))
 
 ngas = made_numerical(gas_FDE, u0)
-T, u, stats = fdae_solver(ngas, [0, 3600], u0.array, Opt(hinit=5, ite_tol=1e-3))
-u = parse_dae_v(u, u0.a)
+T, u, stats = fdae_solver(ngas, [0, 3600], u0, Opt(hinit=5, ite_tol=1e-3))
 ngas_sparse, code = made_numerical(gas_FDE, u0, sparse=True, output_code=True)
 
 T1, u1, stats1 = fdae_solver(ngas_sparse, [0, 3600], u0, Opt(hinit=5, ite_tol=1e-3))
diff --git a/tests/test_gasflow.py b/tests/test_gasflow.py
index 316d810..0ddd69f 100644
--- a/tests/test_gasflow.py
+++ b/tests/test_gasflow.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from Solverz import idx, Var, Para, Sign, as_Vars, nr_method, Eqn, AE, Mat_Mul, made_numerical, parse_ae_v
+from Solverz import idx, Var, Para, Sign, as_Vars, nr_method, Eqn, AE, Mat_Mul, made_numerical
 
 k = idx('k', value=[0, 1, 2])
 i = idx('i', value=[0, 0, 3])
diff --git a/tests/test_m3b9.py b/tests/test_m3b9.py
index e443da5..d8701fb 100644
--- a/tests/test_m3b9.py
+++ b/tests/test_m3b9.py
@@ -2,7 +2,7 @@
 import pandas as pd
 
 from Solverz import DAE, Eqn, Ode, Var, Para, idx, sin, cos, Rodas, as_Vars, Mat_Mul, \
-    implicit_trapezoid, Opt, TimeSeriesParam, made_numerical, parse_dae_v
+    implicit_trapezoid, Opt, TimeSeriesParam, made_numerical
 
 omega = Var(name='omega')
 delta = Var(name='delta')