Repository: BondGraphTools/BondGraphTools (language: python)

Issue:
Cellml edits
Tidying and Bugfixes for CellML Conference presentation


Discussion / hints:

Fix patch:
diff --git a/.flake8 b/.flake8
new file mode 100644
index 0000000..eb31de0
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,3 @@
+[flake8]
+
+ignore = D105, D107
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 74f4a38..8818524 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,5 @@ doctrees/
 build/
 BondGraphTools.egg-info/
 docs/html
+dist
+coverage.xml
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 7297ea5..8cdb3be 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -4,7 +4,7 @@ build:
   image: latest
 
 python:
-  version: 3.6
+  version: 3.7
   setup_py_install: true
   extra_requirements:
     - docs
diff --git a/.travis.yml b/.travis.yml
index 87269d2..0d5a8e2 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,4 +1,4 @@
-dist: xenial
+dist: bionic
 
 env:
   global:
@@ -7,21 +7,23 @@ env:
 language: python
 
 python:
-  - "3.6"
   - "3.7"
+  - "3.8"
+  - "3.9"
 
 before_install:
-  - curl -L https://julialang-s3.julialang.org/bin/linux/x64/0.6/julia-0.6.2-linux-x86_64.tar.gz > julia-0.6.2-linux-x86_64.tar.gz
-  - tar xf julia-0.6.2-linux-x86_64.tar.gz
-  - export PATH=$PATH:./julia-d386e40c17/bin/
-  - chmod +x ./julia-d386e40c17/bin/julia
-
+  - sudo apt-get update && sudo apt-get install -y gfortran libblas-dev liblapack-dev git build-essential cmake python3-dev
+  - git clone https://github.com/LLNL/sundials.git
+  - mkdir sundials/build && cd sundials/build && cmake -DLAPACK_ENABLE=ON ..
+  - make && sudo make install
+  - export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
+  - cd .. && cd ..
 
 install:
+  - echo $PYTHONPATH
+  - pip3 install --upgrade numpy
   - python -m pip install -r requirements.txt
-  - python -m pip install pytest==4.1.0 coverage==4.5.1 mock pytest-cov==2.6.1
-  - LD_LIBRARY_PATH=/usr/local/lib64/:$LD_LIBRARY_PATH
-  - export LD_LIBRARY_PATH
+  - python -m pip install pytest coverage mock pytest-cov==2.6.1
 
 before_script:
   - curl -L https://codeclimate.com/downloads/test-reporter/test-reporter-latest-linux-amd64 > ./cc-test-reporter
diff --git a/BondGraphTools/__init__.py b/BondGraphTools/__init__.py
index f84f558..644843a 100644
--- a/BondGraphTools/__init__.py
+++ b/BondGraphTools/__init__.py
@@ -32,19 +32,19 @@
     resistor = bgt.new("R", name="R1", value=1.0)
 
     # 1 Henry Inductor
-    inductor = bgt.new("L", name="L1", value=1.0)
+    inductor = bgt.new("I", name="L1", value=1.0)
     # 1 Farad Capacitor
     capacitor = bgt.new("C", name="C1", value=1.0)
     # Conservation Law
     law = bgt.new("0") # Common voltage conservation law
 
     # Connect the components
-    connect(law, resistor)
-    connect(law, capacitor)
-    connect(law, inductor)
+    bgt.connect(law, resistor)
+    bgt.connect(law, capacitor)
+    bgt.connect(law, inductor)
 
     # produce timeseries data
-    t, x = simulate(model, x0=[1,1], timespan=[0, 10])
+    t, x = bgt.simulate(model, x0=[1,1], timespan=[0, 10])
 
 Bugs
 ----
diff --git a/BondGraphTools/actions.py b/BondGraphTools/actions.py
index 9357368..bccd371 100644
--- a/BondGraphTools/actions.py
+++ b/BondGraphTools/actions.py
@@ -48,7 +48,7 @@ def disconnect(target, other):
     if isinstance(target, BondGraphBase):
         model = target.parent
     elif isinstance(target, Port):
-        model =  target.component.parent
+        model = target.component.parent
     else:
         model = target[0].parent
 
@@ -100,7 +100,6 @@ def connect(source, destination):
         :func:`disconnect`
     """
 
-
     tail = _find_or_make_port(source, is_tail=True)
     head = _find_or_make_port(destination)
 
@@ -339,8 +338,8 @@ def expose(component, label=None):
             f"Component {component} is not inside anything"
         )
     # fix me with metamodeles or something trickier
-    if component.__component__ is not "SS":
-        ss = new("SS",  name=component.name)
+    if component.__component__ != "SS":
+        ss = new("SS", name=component.name)
         try:
             swap(component, ss)
         except InvalidComponentException as ex:
diff --git a/BondGraphTools/algebra.py b/BondGraphTools/algebra.py
index 2abe43d..7e9b22e 100644
--- a/BondGraphTools/algebra.py
+++ b/BondGraphTools/algebra.py
@@ -18,6 +18,7 @@
     "augmented_rref"
 ]
 
+
 def extract_coefficients(equation: sympy.Expr,
                          local_map: dict,
                          global_coords: list) -> tuple:
@@ -30,13 +31,13 @@ def extract_coefficients(equation: sympy.Expr,
         global_coords: The list of global co-ordinates.
 
     Returns:
-        The linear and nonlinear parts of the equation in the global 
-        co-ordinate system. 
+        The linear and nonlinear parts of the equation in the global
+        co-ordinate system.
 
     Extracts the coordinates from the given equation and maps them into
     the global coordinate space.
     Equations are assumed to come in as sympy expressions of the form
-    :math:`\Phi(x) = 0`.
+    :math:`\\Phi(x) = 0`.
     local_map is a dictionary mappings
 
     .. math::
@@ -93,7 +94,8 @@ def extract_coefficients(equation: sympy.Expr,
     return coeff_dict, nonlinear_terms
 
 
-def _generate_substitutions(linear_op, nonlinear_op, constraints, coords, size_tup):
+def _generate_substitutions(linear_op, nonlinear_op,
+                            constraints, coords, size_tup):
 
     # Lx + F(x) = 0 =>  Ix = (I - L)x - F(x) = Rx - F(x)
     # Since L is in smith normal form (rref and square)
@@ -113,14 +115,14 @@ def _generate_substitutions(linear_op, nonlinear_op, constraints, coords, size_t
     ss_size, js_size, cv_size, n = size_tup
     substitutions = []
     coords_vect = sympy.Matrix(coords)
-    for i in reversed(range(2*(ss_size  + js_size))):
+    for i in reversed(range(2 * (ss_size + js_size))):
         co = coords[i]
-        if Rx[i,i] == 0 and co in atoms and not co in nonlinear_op[i].atoms():
+        if Rx[i, i] == 0 and co in atoms and not co in nonlinear_op[i].atoms():
 
-            eqn = (Rx[i,:]*coords_vect)[0] - nonlinear_op[i]
+            eqn = (Rx[i, :] * coords_vect)[0] - nonlinear_op[i]
             pair = (coords[i], eqn)
             logger.debug("Generating substition %s = %s",
-                        repr(coords[i]), repr(eqn))
+                         repr(coords[i]), repr(eqn))
             substitutions = [
                 (c, s.subs(*pair)) for c, s in substitutions
             ]
@@ -139,20 +141,20 @@ def _process_constraints(linear_op,
     ss_size, js_size, cv_size, n = size_tup
     offset = 2 * js_size + ss_size
 
-    coord_atoms = set(coordinates[0:offset+ss_size])
+    coord_atoms = set(coordinates[0:offset + ss_size])
 
     coords_vect = sympy.Matrix(coordinates)
     cv_constraints = list(
-        linear_op[offset+ss_size:n,:]*coords_vect +
-        nonlinear_op[offset+ss_size:n,0]
+        linear_op[offset + ss_size:n, :] * coords_vect +
+        nonlinear_op[offset + ss_size:n, 0]
     )
     constraints += [cons for cons in cv_constraints if cons]
-    linear_op = linear_op[:offset+ss_size, :]
-    nonlinear_op = nonlinear_op[:offset+ss_size, :]
+    linear_op = linear_op[:offset + ss_size, :]
+    nonlinear_op = nonlinear_op[:offset + ss_size, :]
 
     while constraints:
         constraint, _ = sympy.fraction(constraints.pop())
-        logger.debug("Processing constraint: %s",repr(constraint))
+        logger.debug("Processing constraint: %s", repr(constraint))
         atoms = constraint.atoms() & set(coord_atoms)
 
         # todo: check to see if we can solve f(x) = u => g(u) = x
@@ -168,7 +170,7 @@ def _process_constraints(linear_op,
                     sympy.SparseMatrix(1, linear_op.cols, {(0, idx): 1})
                 )
                 nonlinear_op = nonlinear_op.col_join(
-                    sympy.SparseMatrix(1, 1, {(0,0): -sol})
+                    sympy.SparseMatrix(1, 1, {(0, 0): -sol})
                 )
                 constraint = c - sol
         else:
@@ -178,8 +180,8 @@ def _process_constraints(linear_op,
             partials = [constraint.diff(c) for c in coordinates]
         except Exception as ex:
             logger.exception("Could not differentiate %s with respect to %s",
-                         repr(constraint),repr(coordinates)
-             )
+                             repr(constraint), repr(coordinates)
+                             )
             raise ex
 
         if any(p != 0 for p in partials[0:offset]):
@@ -203,7 +205,7 @@ def _process_constraints(linear_op,
                         nlin += new_coeff * coordinates[idx]
             for idx, coeff in enumerate(cv_derivs):
                 if coeff != 0:
-                    cv = coordinates[offset+ss_size+idx]
+                    cv = coordinates[offset + ss_size + idx]
                     dvc = sympy.Symbol(f"d{str(cv)}")
                     try:
                         dc_idx = coordinates.index(dvc)
@@ -215,16 +217,16 @@ def _process_constraints(linear_op,
                         linear_op = linear_op.row_join(
                             sympy.SparseMatrix(linear_op.rows, 1, {})
                         )
-                    eqn = coeff/factor
+                    eqn = coeff / factor
                     if eqn.is_number:
                         lin_dict.update({(0, dc_idx): eqn})
                     else:
-                        nlin += eqn*dvc
+                        nlin += eqn * dvc
             linear_op = linear_op.col_join(
-                sympy.SparseMatrix(1,linear_op.cols, lin_dict)
+                sympy.SparseMatrix(1, linear_op.cols, lin_dict)
             )
             nonlinear_op = nonlinear_op.col_join(
-                    sympy.SparseMatrix(1,1,{(0,0):nlin})
+                sympy.SparseMatrix(1, 1, {(0, 0): nlin})
             )
 
     linear_op, nonlinear_op, new_constraints = smith_normal_form(
@@ -232,14 +234,14 @@ def _process_constraints(linear_op,
         augment=nonlinear_op)
 
     return linear_op, nonlinear_op, new_constraints + initial_constraints, \
-           coordinates, (ss_size, js_size, cv_size, n)
+        coordinates, (ss_size, js_size, cv_size, n)
 
 
 def _generate_cv_substitutions(subs_pairs, mappins, coords):
     state_map, port_map, control_map = mappins
     ss_size = len(state_map)
 
-    cv_offset = 2*(ss_size + len(port_map))
+    cv_offset = 2 * (ss_size + len(port_map))
 
     control_vars = {str(c) for c in coords[cv_offset:]}
     subs = []
@@ -341,30 +343,32 @@ def reduce_model(linear_op, nonlinear_op, coordinates, size_tuple,
             for idx in range(ss_size):
                 v = linear_op[row, idx + offset]
                 if v:
-                    lin_dict.update({(rows_added,idx): v})
+                    lin_dict.update({(rows_added, idx): v})
             for idx in range(cv_size):
-                v = linear_op[row, idx + offset+ss_size]
+                v = linear_op[row, idx + offset + ss_size]
                 if v:
                     cv_diff_dict.update({(rows_added, idx): v})
 
     for row in range(offset, linear_op.rows):
         logger.debug("Testing row %s: %s + %s", repr(row),
-                    repr(linear_op[row, :] * sympy.Matrix(coordinates)),
-                    repr(nonlinear_op[row]) if nonlinear_op else '')
+                     repr(linear_op[row, :] * sympy.Matrix(coordinates)),
+                     repr(nonlinear_op[row]) if nonlinear_op else '')
 
         nonlinear_constraint = nonlinear_op[row]
         F_args = set(coordinates[0:offset + ss_size]) & \
-                 nonlinear_constraint.atoms()
-        if linear_op[row, offset:-1].is_zero and not nonlinear_constraint:
+            nonlinear_constraint.atoms()
+        if linear_op[row, offset:-
+                     1].is_zero_matrix and not nonlinear_constraint:
             continue
 
         state_constraint = linear_op[row, offset: offset + ss_size]
         control_constraint = linear_op[row, offset + ss_size:]
 
-        row = state_constraint.row_join(sympy.SparseMatrix(1, offset + cv_size, {}))
+        row = state_constraint.row_join(
+            sympy.SparseMatrix(1, offset + cv_size, {}))
 
         cv_dict = {}
-        if not control_constraint.is_zero:
+        if not control_constraint.is_zero_matrix:
             logger.debug("Found higher order control constraint")
             for cv_col in range(control_constraint.cols):
                 const = control_constraint[cv_col]
@@ -376,14 +380,15 @@ def reduce_model(linear_op, nonlinear_op, coordinates, size_tuple,
                 except ValueError:
                     idx = len(added_cvs)
                     added_cvs.append(cv_col)
-                    linear_op= linear_op.row_join(sympy.SparseMatrix(linear_op.rows, 1, {}))
+                    linear_op = linear_op.row_join(
+                        sympy.SparseMatrix(linear_op.rows, 1, {}))
                     coord = coordinates[offset + ss_size + cv_col]
                     d_coord = sympy.Symbol(f"d{str(coord)}")
                     coordinates.append(d_coord)
                     cv_size += 1
                     n += 1
 
-                cv_dict[(0,idx)] = const
+                cv_dict[(0, idx)] = const
 
         row = row.row_join(sympy.SparseMatrix(1, len(added_cvs), cv_dict))
         jac_dx = [nonlinear_constraint.diff(c) for c in coordinates[:ss_size]]
@@ -394,7 +399,7 @@ def reduce_model(linear_op, nonlinear_op, coordinates, size_tuple,
         jac_x = [
             nonlinear_constraint.diff(c)
             for c in coordinates[offset:
-            offset+ss_size]
+                                 offset + ss_size]
         ]
         jac_cv = [
             nonlinear_constraint.diff(c)
@@ -403,26 +408,26 @@ def reduce_model(linear_op, nonlinear_op, coordinates, size_tuple,
 
         nlin_row = sympy.S(0)
 
-        if any(x!=0 for x in jac_dx):
+        if any(x != 0 for x in jac_dx):
             logger.warning("Second order constraint not implemented: %s",
                            jac_dx)
 
-        elif any(x!=0 for x in jac_junciton):
+        elif any(x != 0 for x in jac_junciton):
             logger.warning("First order junciton constraint not implemented: %s",
                            str(jac_junciton))
 
-        elif any(x!=0 for x in jac_cv):
+        elif any(x != 0 for x in jac_cv):
             logger.warning("First order control constraint not implemented: %s",
                            str(jac_cv))
 
-        elif any(x!=0 for x in jac_x):
+        elif any(x != 0 for x in jac_x):
             logger.debug("First order constriants: %s", jac_x)
-            fx = sum(x*y for x,y in zip(jac_x, coordinates[:ss_size]))
+            fx = sum(x * y for x, y in zip(jac_x, coordinates[:ss_size]))
             logger.debug(repr(fx))
             p, q = sympy.fraction(sympy.simplify(fx))
-            if row.is_zero:
+            if row.is_zero_matrix:
                 lin_dict, nlin = extract_coefficients(
-                    p, {c:i for i,c in enumerate(coordinates)},
+                    p, {c: i for i, c in enumerate(coordinates)},
                     coordinates)
 
                 for k, v in lin_dict.items():
@@ -433,7 +438,8 @@ def reduce_model(linear_op, nonlinear_op, coordinates, size_tuple,
             else:
                 nlin_row += fx
 
-        nonlinear_op = nonlinear_op.col_join(sympy.SparseMatrix(1,1,[nlin_row]))
+        nonlinear_op = nonlinear_op.col_join(
+            sympy.SparseMatrix(1, 1, [nlin_row]))
 
         linear_op = linear_op.col_join(row)
         rows_added += 1
@@ -489,7 +495,7 @@ def augmented_rref(matrix, augmented_rows=0):
                     if v > v_max:
                         j = row
                         v_max = v
-                except TypeError: # symbolic variable
+                except TypeError:  # symbolic variable
                     j = row
                     v_max = v
             if not j:
@@ -501,7 +507,7 @@ def augmented_rref(matrix, augmented_rows=0):
 
         for i in range(matrix.rows):
             if i != pivot and matrix[i, col] != 0:
-                b = matrix[i, col]/a
+                b = matrix[i, col] / a
                 matrix[i, :] += - b * matrix[pivot, :]
 
         matrix[pivot, :] *= 1 / a
@@ -558,25 +564,25 @@ def smith_normal_form(matrix, augment=None):
     m, n = M.shape
     M = augmented_rref(M, k)
 
-    Mp = sympy.MutableSparseMatrix(n-k, n, {})
+    Mp = sympy.MutableSparseMatrix(n - k, n, {})
 
     constraints = []
     for row in range(m):
         leading_coeff = -1
-        for col in range(row, n-k):
+        for col in range(row, n - k):
             if M[row, col] != 0:
                 leading_coeff = col
                 break
         if leading_coeff < 0:
-            if not M[row, n-k:].is_zero:
-                constraints.append(sum(M[row,:]))
+            if not M[row, n - k:].is_zero_matrix:
+                constraints.append(sum(M[row, :]))
         else:
             Mp[leading_coeff, :] = M[row, :]
 
     if augment:
-        return Mp[:,:-k], Mp[:, -k:], constraints
+        return Mp[:, :-k], Mp[:, -k:], constraints
     else:
-        return Mp, sympy.SparseMatrix(m,k,{}), constraints
+        return Mp, sympy.SparseMatrix(m, k, {}), constraints
 
 
 def adjacency_to_dict(nodes, edges, offset=0):
@@ -643,20 +649,20 @@ def get_relations_iterator(component, mappings, coordinates, io_map=None):
     # todo: Fix this dirty hack; there has to be a better way to hand io ports
     for cv, value in local_cv.items():
         try:
-            local_map[cv] = 2*(num_ports+num_state_vars) + inv_cv[value]
+            local_map[cv] = 2 * (num_ports + num_state_vars) + inv_cv[value]
         except KeyError:
             logger.debug("Could not find %s, trying the io_map", value)
             key = io_map[value]
             local_map[cv] = key
-            logger.debug("Mapping %s to co-ord %s",cv, coordinates[key])
+            logger.debug("Mapping %s to co-ord %s", cv, coordinates[key])
 
     for (x, dx), coord in local_tm.items():
         local_map[dx] = inv_tm[coord]
         local_map[x] = inv_tm[coord] + num_state_vars + 2 * num_ports
 
     for (e, f), port in local_js.items():
-        local_map[e] = 2*inv_js[port] + num_state_vars
-        local_map[f] = 2*inv_js[port] + num_state_vars + 1
+        local_map[e] = 2 * inv_js[port] + num_state_vars
+        local_map[f] = 2 * inv_js[port] + num_state_vars + 1
     logger.debug("Getting relations iterator for %s", repr(component))
     for relation in component.constitutive_relations:
         if relation:
diff --git a/BondGraphTools/atomic.py b/BondGraphTools/atomic.py
index e0e0ada..4ca2ecb 100644
--- a/BondGraphTools/atomic.py
+++ b/BondGraphTools/atomic.py
@@ -13,9 +13,7 @@
 
 
 class Component(BondGraphBase, PortManager):
-    """
-    Atomic bond graph components are those defined by constitutive relations.
-    """
+    """Components defined by constitutive relations."""
 
     def __init__(self, metamodel, constitutive_relations,
                  state_vars=None, params=None, **kwargs):
@@ -180,8 +178,8 @@ def _build_relations(self):
                 for port_id in self.ports:
                     if isinstance(port_id, int):
                         rels.append(
-                        sp.sympify(
-                            string.replace("_i", "_{}".format(port_id))))
+                            sp.sympify(
+                                string.replace("_i", "_{}".format(port_id))))
             else:
 
                 tiers = 0
@@ -203,7 +201,7 @@ def _build_relations(self):
                 terms = [substring.replace("_i", "_{}".format(p))
                          for p in self.ports if isinstance(p, int)]
                 symstr = string[0:sloc] + "(" + " + ".join(terms) + string[
-                                                                    eloc:]
+                    eloc:]
 
                 rels.append(
                     sp.sympify(symstr)
@@ -222,6 +220,7 @@ class SymmetricComponent(Component):
     Instances of this class are multi-port components which are able to have
     connections made without specifying ports.
     """
+
     def get_port(self, port=None):
         """See `Component`"""
         if not port and not isinstance(port, int):
@@ -245,6 +244,7 @@ class EqualEffort(BondGraphBase, PortExpander):
          constitutive_relations:
 
     """
+
     def __init__(self, **kwargs):
 
         PortExpander .__init__(self, {None: None})
@@ -304,7 +304,7 @@ def template(self):
 
     @property
     def basis_vectors(self):
-        return {},  self._port_vectors(), {}
+        return {}, self._port_vectors(), {}
 
     def get_port(self, port=None):
         try:
@@ -322,13 +322,13 @@ def constitutive_relations(self):
         (e_0, f_0), port = var.pop()
 
         sigma_0 = port.weight
-        partial_sum = sigma_0*e_0
+        partial_sum = sigma_0 * e_0
 
         while var:
             (e_i, f_i), port = var.pop()
             sigma_i = port.weight
-            partial_sum += sigma_i*e_i
-            relations.append(sigma_i*f_i - sigma_0*f_0)
+            partial_sum += sigma_i * e_i
+            relations.append(sigma_i * f_i - sigma_0 * f_0)
 
         relations.append(partial_sum)
         return relations
diff --git a/BondGraphTools/base.py b/BondGraphTools/base.py
index 1d5f388..46b355b 100644
--- a/BondGraphTools/base.py
+++ b/BondGraphTools/base.py
@@ -1,5 +1,5 @@
-"""This module contains the base classes for bond graph models and connections
-"""
+"""Base classes for bond graph models and connections."""
+
 import logging
 from collections import namedtuple
 
@@ -10,19 +10,11 @@ class BondGraphBase(object):
     """
     Base class definition for all bond graphs.
 
-    Attributes:
-        parent:
-        name:
-        metamodel:
-        template:
-        uri:
-        root:
-        basis_vectors:
-
     Args:
-        name: Assumed to be unique
-        parent:
-        metadata (dict):
+        name:               Name of this model, assumed to be unique.
+        parent:             Parent model, or none.
+        metamodel:          Metamodel class.
+
     """
 
     def __init__(self,
@@ -44,12 +36,13 @@ def __init__(self,
             self.name = name
 
         self.parent = parent
+        """Model that contains this object."""
         self.view = None
 
     def __repr__(self):
         return f"{self.metamodel}: {self.name}"
 
-    def __new__(cls, *args, **kwargs):
+    def __new__(cls, *args, **kwargs):  # noqa: D102
         if "instances" not in cls.__dict__:
             cls.instances = 1
         else:
@@ -62,18 +55,22 @@ def __del__(self):
 
     @property
     def metamodel(self):
+        """The meta-model type of this object."""  # noqa: D401
         return self.__metamodel
 
     @property
     def template(self):
+        """The model template from which this was created."""  # noqa: D401
         raise NotImplementedError
 
     @property
     def constitutive_relations(self):
+        """The equations governing this objects dynamics."""   # noqa: D401
         raise NotImplementedError
 
     @property
     def uri(self):
+        """Model reference locator."""  # noqa: D401
         if not self.parent:
             return f"{self.name}:"
         else:
@@ -81,6 +78,7 @@ def uri(self):
 
     @property
     def root(self):
+        """The root of the tree to which this object belongs."""  # noqa: D401
         if not self.parent:
             return self
         else:
@@ -88,6 +86,7 @@ def root(self):
 
     @property
     def basis_vectors(self):
+        """The input/ouput, dynamic and interconnect vectors."""  # noqa: D401
         raise NotImplementedError
 
     def __hash__(self):
@@ -95,13 +94,16 @@ def __hash__(self):
 
 
 class Bond(namedtuple("Bond", ["tail", "head"])):
-    """A `namedtuple` that stores a connection between two ports.
+    """Stores the connection between two ports.
+
     Head and tail are specified to determine orientation
 
     Attributes:
         head: The 'harpoon' end of the power bond and direction of positive $f$
         tail: The non-harpoon end, and direction of negative $f$
+
     """
+
     def __contains__(self, item):
         if isinstance(item, BondGraphBase):
             return self.head[0] is item or self.tail[0] is item
@@ -114,17 +116,22 @@ def __contains__(self, item):
 
 
 class Port(object):
-    """
-    Basic object for representing ports;
+    """Basic object for representing ports.
+
     Looks and behaves like a namedtuple:
-    component, index =  Port
+        component, index =  Port
+
+    Args:
+        component:      The owner of this port.
+        index:          The index of this port in the component.
+
     """
 
     def __init__(self, component, index):
         self.component = component
         """(`PortManager`) The component that this port is attached to"""
         self.index = index
-        """(int) The numberical index of this port"""
+        """(int) The numerical index of this port"""
         self.is_connected = False
         """(bool) True if this port is plugged in."""
 
@@ -159,7 +166,7 @@ def __repr__(self):
 
     def __eq__(self, other):
         try:
-            return ((self.component is other.component) and
+            return ((self.component is other.component) and  # noqa
                     (self.index == other.index))
         except AttributeError:
             try:
diff --git a/BondGraphTools/component_manager.py b/BondGraphTools/component_manager.py
index 7ac55e4..67a3091 100644
--- a/BondGraphTools/component_manager.py
+++ b/BondGraphTools/component_manager.py
@@ -1,4 +1,5 @@
-"""
+"""Component Libraries.
+
 This module takes care of the loading and management of component libraries.
 This library will automatically load all factory libraries upon import.
 Additionally libraries can be added via :func:`load_library`.
@@ -36,14 +37,14 @@
 
 
 def load_library(filename):
-    """
-    Loads the library specified by the filename and makes it available for use.
+    """Load the library specified by the filename.
 
     Args:
         filename: The (absolute) filename of the library to be loaded.
 
     Returns:
         bool: True if the library was successfully loaded.
+
     """
     with open(filename) as file:
         try:
@@ -69,34 +70,34 @@ def load_library(filename):
 
 
 def get_library_list():
-    """
-    Fetches a list of the libraries available for use.
+    """Fetch a list of the libraries available for use.
 
     Returns:
         list of (library id, description) tuples
-    """
 
-    return [(l, __libraries[l]["description"]) for l in __libraries]
+    """
+    return [(lib, __libraries[lib]["description"]) for lib in __libraries]
 
 
 def get_components_list(library):
-    """
-    Fetches a list of components available in the given library.
+    """Fetch a list of components available in the given library.
 
     Args:
         library: The library id of the library to query
 
     Returns:
         list of (component id, description) tuples
+
     """
     components = __libraries[library]["components"]
 
-    return [(comp_id, components[comp_id]["description"]) for comp_id in components]
+    return [(comp_id, components[comp_id]["description"])
+            for comp_id in components]
 
 
 def get_component(component, library=base_id):
-    """
-    Fetches the component data for the specified component
+    """Fetch the component data for the specified component.
+
     Args:
         component: The id of the specific component
         library: The id of the library to which the component belongs
@@ -109,42 +110,40 @@ def get_component(component, library=base_id):
 
 
 def find(component, restrict_to=None, find_all=False, ensure_unique=False):
-    """
-    Finds the specified component.
+    """Find the specified component.
 
     Args:
-        component: The component id to find.
-        restrict_to: `list` or `set` of library id's to be that the search
-         should be restricted to.
-        find_all: `False` if the function should return only the first instance
-         of the component, `True` if the function should return all such
-         instances
-        ensure_unique: If true, this assumes that the component id must be
-         unique across libraries, and hence will raise an exception if this is
-         assumption is violated.
-
-    Returns: the library id, or a list of library_id in which this component
-     can be found.
+        component:     The component id to find.
+        restrict_to:   `list` or `set` of library id's to be that the search
+                       should be restricted to.
+        find_all:      `False` if the function should return only the first
+                        instance of the component, `True` if the function
+                        should return all such instances
+        ensure_unique:  If true, this assumes that the component id must be
+                        unique across libraries, and hence will raise an
+                        exception if this is assumption is violated.
+
+    Returns:
+        the library id, or a list of library_id in which this component
+        can be found.
 
     Raises:
         NotImplementedError - if the component is not found.
         ValueError - if the component is assume to be unique, but is not
-    """
-    results = []
 
+    """
     if restrict_to:
         keys = {k for k in __libraries if k in restrict_to}
     else:
         keys = __libraries
 
-    unique_id = None
-
-    results = {l for l in keys if component in __libraries[l]["components"]}
+    results = {lib for lib in keys
+               if component in __libraries[lib]["components"]}
 
     if find_all:
         return results
     elif (ensure_unique and len(results) == 1) or not ensure_unique:
-        return results[0]
+        return list(results)[0]
     elif ensure_unique:
         raise KeyError("Component is not unique")
     else:
diff --git a/BondGraphTools/compound.py b/BondGraphTools/compound.py
index 97c8ed2..5c25da5 100644
--- a/BondGraphTools/compound.py
+++ b/BondGraphTools/compound.py
@@ -23,12 +23,13 @@
 class BondGraph(BondGraphBase, LabeledPortManager):
     """Representation of a bond graph model.
     """
+
     def __init__(self, name, components=None, **kwargs):
 
         BondGraphBase.__init__(self, name, **kwargs)
         LabeledPortManager.__init__(self)
         self.components = OrderedSet()
-        """The components, instances of :obj:`BondGraphBase`, 
+        """The components, instances of :obj:`BondGraphBase`,
         that make up this model"""
 
         if components:
@@ -225,7 +226,7 @@ def control_vars(self):
         out = dict()
         excluded = {
             v for pair in self._port_map.values() for v in pair
-         }
+        }
 
         for v in self.components:
             try:
@@ -302,8 +303,8 @@ def constitutive_relations(self):
         coord_vect = sp.Matrix(coordinates)
         relations = [
             sp.Add(l, r) for i, (l, r) in enumerate(zip(
-                lin_op*coord_vect, nlin_op))
-            if not ss_size <= i < ss_size + 2*js_size - 2*len(out_ports)
+                lin_op * coord_vect, nlin_op))
+            if not ss_size <= i < ss_size + 2 * js_size - 2 * len(out_ports)
         ]
         if isinstance(constraints, list):
             for constraint in constraints:
@@ -344,8 +345,8 @@ def system_model(self, control_vars=None):
 
         .. math::
 
-            x = (dx_0,  dx_1,  \ldots,  e_0,  f_0, e_1, f_1, \ldots, x_0,
-            x_1, \ldots, u_0, u_1, \ldots)
+            x = (dx_0,  dx_1,  \\ldots,  e_0,  f_0, e_1, f_1, \\ldots, x_0,
+            x_1, \\ldots, u_0, u_1, \\ldots)
 
         So that the system obeys the differential-algebraic equation
 
@@ -382,8 +383,8 @@ def system_model(self, control_vars=None):
         inverse_port_map = {}
 
         for port, (cv_e, cv_f) in self._port_map.items():
-            inverse_port_map[cv_e] = ss_size + 2*inv_js[port]
-            inverse_port_map[cv_f] = ss_size + 2*inv_js[port] + 1
+            inverse_port_map[cv_e] = ss_size + 2 * inv_js[port]
+            inverse_port_map[cv_f] = ss_size + 2 * inv_js[port] + 1
 
         for component in self.components:
             relations = get_relations_iterator(
@@ -399,7 +400,7 @@ def system_model(self, control_vars=None):
         linear_op = sp.SparseMatrix(row, n, lin_dict)
         nonlinear_op = sp.SparseMatrix(row, 1, nlin_dict)
         coordinates, linear_op, nonlinear_op, constraints = reduce_model(
-                linear_op, nonlinear_op, coordinates, size_tuple,
+            linear_op, nonlinear_op, coordinates, size_tuple,
             control_vars=control_vars
         )
 
@@ -459,6 +460,7 @@ class BondSet(OrderedSet):
     """
     Container class for internal bonds.
     """
+
     def add(self, bond):
         tail = bond.tail
         head = bond.head
diff --git a/BondGraphTools/config.py b/BondGraphTools/config.py
deleted file mode 100644
index 68828ce..0000000
--- a/BondGraphTools/config.py
+++ /dev/null
@@ -1,180 +0,0 @@
-"""Manages system configuration
-
-This file handles the system configuration, and in particular the bridge between
-python and julia.
-
-To access the julia interpreter use::
-    julia = config.julia
-
-"""
-
-import pathlib
-import json
-import sys
-import os
-import logging
-from subprocess import Popen, PIPE, run
-
-from .version import __version__ as VERSION
-logger = logging.getLogger(__name__)
-
-
-def _as_str(path):
-    if isinstance(path, str):
-        return path
-    else:
-        return str(path.as_posix())
-
-
-def _check_julia(command):
-    # Check to see if Julia is already in the path
-
-    p = Popen([command, "-v"], stdout=PIPE, stderr=None)
-
-    out, err = p.communicate()
-    text = out.decode('utf8').strip()
-    if text.startswith("julia version 0.6"):
-        return True
-    else:
-        return False
-
-
-class Config:
-    _which = 'which'
-    base = pathlib.Path.home().absolute() / '.BondGraphTools'
-    file = base / 'config.json'
-
-    def __init__(self, julia_executable=None, python_executable=None, **kwargs):
-
-        self._julia = None
-        self._de = None
-        new_config = False
-        rebuild = False
-
-        if not julia_executable:
-            self.julia_executable = self.find_julia()
-            new_config = True
-        else:
-            self.julia_executable = pathlib.Path(julia_executable)
-
-        if python_executable and python_executable is not sys.executable:
-            rebuild = True
-        elif python_executable:
-            self.python_executable = pathlib.Path(python_executable)
-        else:
-            self.python_executable = pathlib.Path(sys.executable).resolve()
-
-        if new_config:
-            if not self.base.exists():
-                self.base.mkdir()
-                self.file.touch()
-            self.install_dependencies(rebuild)
-            self.save()
-
-    def find_julia(self):
-        if _check_julia('julia'):
-            p = Popen([self._which, 'julia'], stdout=PIPE, stderr=None)
-            out, err = p.communicate()
-            julia_dir = out.decode('utf8').strip()
-            return pathlib.Path(julia_dir).resolve()
-        else:
-            raise NotImplementedError("Could not find Julia 0.6.4: "
-                                      "please install and add to path")
-
-    def find_conda(self):
-        p = Popen([self._which, 'conda'], stdout=PIPE, stderr=None)
-        out, err = p.communicate()
-        conda_dir = out.decode('utf8').strip()
-        if not conda_dir:
-            return None
-        else:
-            return pathlib.Path(conda_dir).resolve()
-
-    def install_dependencies(self, rebuild=False):
-        # we assume julia and python are already in the path
-        logger.warning('Installing Julia dependencies; '
-                       'this may take some time')
-
-        env = os.environ
-        env.update({"PYTHON": _as_str(self.python_executable),
-                    "JULIA": _as_str(self.julia_executable)})
-        conda = self.find_conda()
-        if conda:
-            env.update({"CONDA": _as_str(conda)})
-
-        julia_code = [
-            "Pkg.init()\n",
-            "Pkg.update()\n",
-            """Pkg.add("PyCall")\n""",
-            """Pkg.add("DifferentialEquations")\n"""]
-
-        if rebuild:
-            julia_code += [
-                """Pkg.build("PyCall")""",
-                "using PyCall\n",
-                "using DifferentialEquations\n"]
-
-        temp = self.base / 'deps.jl'
-        temp.touch()
-        with open(_as_str(temp), 'w') as fs:
-            fs.writelines(julia_code)
-        run([_as_str(self.julia_executable), _as_str(temp)], env=env)
-        os.remove(temp)
-        logger.warning("Complete")
-
-    @property
-    def julia(self):
-        if not self._julia:
-            self.start_julia()
-
-        return self._julia
-
-    def start_julia(self):
-
-        import diffeqpy
-        from diffeqpy import de
-
-        self._julia = diffeqpy.setup()
-        self._de = de
-
-    @property
-    def de(self):
-        if not self._julia:
-            self.start_julia()
-
-        return self._de
-
-    @staticmethod
-    def load():
-        try:
-            with open(_as_str(Config.file), 'r') as f:
-                kwargs = json.load(f)
-        except (json.JSONDecodeError, FileNotFoundError) as ex:
-            kwargs = {}
-
-        if sys.platform.startswith('win'):
-            return WinConfig(**kwargs)
-        else:
-            return Config(**kwargs)
-
-    def save(self):
-
-        config_dict = {
-            'julia_executable': _as_str(self.julia_executable),
-            'python_executable': _as_str(self.python_executable),
-            'version': VERSION
-        }
-        with open(_as_str(Config.file), 'w') as f:
-            json.dump(config_dict, f)
-
-
-class WinConfig(Config):
-    """
-    Extra configuration functionality for windows setup.
-    """
-    _which = 'where'
-
-
-config = Config.load()
-
-__all__ = ["config"]
diff --git a/BondGraphTools/fileio.py b/BondGraphTools/fileio.py
index 551c4a9..2e56174 100644
--- a/BondGraphTools/fileio.py
+++ b/BondGraphTools/fileio.py
@@ -1,4 +1,4 @@
-"""The file save/load interface and file format data model
+"""The file save/load interface and file format data model.
 
 This module provides the basic IO functionality such as saving and loading
 to file.
@@ -14,21 +14,21 @@
 
 from .compound import BondGraph
 from .actions import connect, new, expose
-from .exceptions import *
+from .exceptions import InvalidComponentException
 
 logger = logging.getLogger(__name__)
 
 FILE_VERSION = "0.1"
 
+
 def save(model, filename):
-    """
-    Saves the model to the specified path
+    """Save the model to file.
 
     Args:
         model: The model to be saved
         filename: The file to save to
-    """
 
+    """
     model_directory = _build_model_directory(model)
 
     models = {}
@@ -39,7 +39,7 @@ def save(model, filename):
             uri: _build_model_data(sub_model, templates)
         })
     data = {
-        "version":FILE_VERSION,
+        "version": FILE_VERSION,
         "root": model.name,
         "models": models
     }
@@ -89,8 +89,8 @@ def _build_model_data(model, templates):
 
     ports = []
     for port in model.ports:
-        (c1, e), (c2, f) = model._port_map[port]
-        if c1==c2:
+        (c1, e), (c2, f) = model._port_map[port]  # noqa
+        if c1 == c2:
             ports.append(f"{c1.name} {port.name}")
         else:
             raise NotImplementedError
@@ -120,13 +120,12 @@ def _build_component_string(component):
 
     except AttributeError:
         pass
-    logger.debug("Saving component string: %s", out_str )
+    logger.debug("Saving component string: %s", out_str)
     return out_str
 
 
 def load(file_name, model=None, as_name=None):
-    """
-    Loads a model from file
+    """Load a model from file.
 
     Args:
         file_name (str or Path): The file to load.
@@ -136,12 +135,13 @@ def load(file_name, model=None, as_name=None):
 
     Raises:
         `NotImplementedError` for incorrect file version.
+
     """
     if isinstance(file_name, pathlib.Path):
         file_name = str(file_name)
 
     with open(file_name, 'r') as f:
-        data = yaml.load(f)
+        data = yaml.load(f, Loader=yaml.SafeLoader)
 
     version = str(data['version'])
     if version == FILE_VERSION:
@@ -180,8 +180,8 @@ def _build(model_name, template_name):
         _wire(model, netlist)
 
         try:
-            IO_ports = model_data["ports"]
-            _expose(model, IO_ports)
+            io_ports = model_data["ports"]
+            _expose(model, io_ports)
         except KeyError:
             logger.debug("No ports on model ")
 
@@ -205,7 +205,7 @@ def get_port(port_string):
             except StopIteration:
                 return comp
             try:
-                t3 = next(tokens)
+                _ = next(tokens)
             except StopIteration:
 
                 logger.debug("Tyring to get port %s, %s", str(comp), str(t2))
@@ -227,8 +227,8 @@ def get_port(port_string):
             head = get_port(head_str)
             connect(tail, head)
 
-    def _expose(model, IO_ports):
-        for port_string in IO_ports:
+    def _expose(model, _io_ports):
+        for port_string in _io_ports:
             component_name, port_label = port_string.split(" ")
             comp, = {
                 c for c in model.components if c.name == component_name
@@ -238,7 +238,7 @@ def _expose(model, IO_ports):
     def _parse_build_args(in_args):
 
         if not in_args:
-            return [] , {}
+            return [], {}
 
         arg = in_args[0]
         args, kwargs = _parse_build_args(in_args[1:])
@@ -253,11 +253,11 @@ def _parse_build_args(in_args):
             try:
                 v = int(v)
             except ValueError:
-                v= float(v)
+                v = float(v)
         if not k:
             args.append(v)
         else:
-            kwargs.update({k:v})
+            kwargs.update({k: v})
 
         return args, kwargs
 
@@ -275,13 +275,6 @@ def _base_component_from_str(string):
             comp.set_param(k, v)
         return comp
 
-    def _validate(model_dict):
-        # We should use this function to make sure we don't have an infinite loop
-        # in the build cycle.
-        pass
-
-    _validate(models)
-
     out = _build(root, root)
 
     if as_name:
@@ -290,6 +283,3 @@ def _validate(model_dict):
         out.name = data["root"]
 
     return out
-
-
-
diff --git a/BondGraphTools/port_hamiltonian.py b/BondGraphTools/port_hamiltonian.py
index 83b99a0..a4526e9 100644
--- a/BondGraphTools/port_hamiltonian.py
+++ b/BondGraphTools/port_hamiltonian.py
@@ -29,7 +29,7 @@ class PortHamiltonian(Component):
 
         creates a new port hamiltonian component `ph`, which has one port
         defined by the relation
-        $$\dot{x}_i = f_i, \qquad e_i = \frac{\partial H}{\partial x_i}$$
+        $$\\dot{x}_i = f_i, \\qquad e_i = \frac{\\partial H}{\\partial x_i}$$
         where $H = w*x^2$,.
 
     See Also: Component
diff --git a/BondGraphTools/port_managers.py b/BondGraphTools/port_managers.py
index 15584a1..b378bb1 100644
--- a/BondGraphTools/port_managers.py
+++ b/BondGraphTools/port_managers.py
@@ -92,6 +92,7 @@ class PortExpander(PortManager):
 
     See Also: `PortTemplate`
     """
+
     def __init__(self, ports, static_ports=None):
         if static_ports:
             super().__init__(static_ports)
@@ -211,6 +212,7 @@ def spawn(self, index=None):
 
 class LabeledPort(Port):
     """See Also: `Port`, `LabeledPortManager`"""
+
     def __init__(self, *args, name=None, **kwargs):
         self.name = name
         """The name of this port"""
@@ -232,6 +234,7 @@ class LabeledPortManager(PortManager):
     See Also: `PortManager`
 
     """
+
     def __init__(self, ports=None):
         if ports:
             super().__init__(ports)
diff --git a/BondGraphTools/reaction_builder.py b/BondGraphTools/reaction_builder.py
index 9f188c5..769ba27 100644
--- a/BondGraphTools/reaction_builder.py
+++ b/BondGraphTools/reaction_builder.py
@@ -28,6 +28,7 @@ class Reaction_Network(object):
         temperature: Temperature in Kelvin (Default 300K, or approx 27c)
         volume:
     """
+
     def __init__(self, reactions=None, name=None, temperature=300, volume=1):
 
         self._reactions = {}
@@ -60,7 +61,7 @@ def stoichiometry(self):
     @property
     def reverse_stoichiometry(self):
         """The reverse stoichiometric matrix"""
-        matrix = SparseMatrix(len(self._species), len(self._reactions),{})
+        matrix = SparseMatrix(len(self._species), len(self._reactions), {})
         for col, (_, forward_species, _, _) in enumerate(
                 self._reactions.values()):
             for species, qty in forward_species.items():
@@ -71,8 +72,8 @@ def reverse_stoichiometry(self):
     @property
     def forward_stoichiometry(self):
         """The forward stoichiometric matrix"""
-        matrix = SparseMatrix(len(self._species), len(self._reactions),{})
-        for col, (back_species,_, _, _) in enumerate(
+        matrix = SparseMatrix(len(self._species), len(self._reactions), {})
+        for col, (back_species, _, _, _) in enumerate(
                 self._reactions.values()):
             for species, qty in back_species.items():
                 matrix[(self.species.index(species), col)] = qty
@@ -108,10 +109,10 @@ def fluxes(self):
                             symbols(f"x_{{{s}}}"), qty)
             eqn = prod - subst
             fluxes.append(eqn)
-        matrix = Matrix(len(fluxes),1, fluxes)
+        matrix = Matrix(len(fluxes), 1, fluxes)
         return matrix, coords
 
-    def as_network_model(self, normalised: bool=False):
+    def as_network_model(self, normalised: bool = False):
         """Produces a bond graph :obj:`BondGraph.BondGraph` model of the system
 
         Args:
@@ -175,18 +176,18 @@ def _connect_complex(self, system, species_anchors, junct, stoichiometry,
 
             if qty == 1:
                 if is_reversed:
-                    connect(species_anchors[species],junct.non_inverting)
+                    connect(species_anchors[species], junct.non_inverting)
                 else:
-                    connect(junct.non_inverting,species_anchors[species])
+                    connect(junct.non_inverting, species_anchors[species])
             else:
                 tf = new("TF", value=qty)
                 system.add(tf)
                 if is_reversed:
                     connect((tf, 1), junct.non_inverting)
-                    connect(species_anchors[species], (tf,0))
+                    connect(species_anchors[species], (tf, 0))
                 else:
                     connect(junct.non_inverting, (tf, 1))
-                    connect((tf,0), species_anchors[species])
+                    connect((tf, 0), species_anchors[species])
 
     def _build_species(self, system, normalised):
         if normalised:
@@ -199,12 +200,12 @@ def _build_species(self, system, normalised):
             # edit: create new component for chemostat
             if species in self._chemostats:
                 this_species = new(
-                        "Se", name=species, value=self._chemostats[species]
+                    "Se", name=species, value=self._chemostats[species]
                 )
-                n_reactions = n_reactions-1
+                n_reactions = n_reactions - 1
             else:
                 this_species = new(
-                        "Ce", library=LIBRARY, name=species, value=param_dict
+                    "Ce", library=LIBRARY, name=species, value=param_dict
                 )
             system.add(this_species)
 
@@ -260,7 +261,7 @@ def add_reaction(self, reaction,
 
         if not name or name in self._reactions:
             n = 1
-            while "r_{{{base}}}".format(base=str(n)+';0') in self._reactions:
+            while "r_{{{base}}}".format(base=str(n) + ';0') in self._reactions:
                 n += 1
             base = str(n) + ';'
             idx = "r_{{{base}}}"
@@ -288,12 +289,12 @@ def add_reaction(self, reaction,
             for sp in reaction_step[i]:
                 self._species[sp] += 1
 
-            for sp in reaction_step[i+1]:
+            for sp in reaction_step[i + 1]:
                 self._species[sp] += 1
 
-            self._reactions[idx.format(base=base+str(i))] = (reaction_step[i],
-                                              reaction_step[i+1],
-                                              f_rate, r_rate)
+            self._reactions[idx.format(base=base + str(i))] = (reaction_step[i],
+                                                               reaction_step[i + 1],
+                                                               f_rate, r_rate)
 
     def add_chemostat(self, species, concentration=None):
         """
@@ -355,4 +356,3 @@ def _split_reactants(reactants):
         stoiciometrics[prod] = coeff
 
     return stoiciometrics
-
diff --git a/BondGraphTools/sim_tools.py b/BondGraphTools/sim_tools.py
index 0ce44f7..80c5a2d 100644
--- a/BondGraphTools/sim_tools.py
+++ b/BondGraphTools/sim_tools.py
@@ -7,13 +7,13 @@
 import sympy as sp
 from sympy.core import SympifyError
 from scipy.optimize import broyden1
-
+from scikits.odes.dae import dae
 from .exceptions import ModelException, SolverException
 
 logger = logging.getLogger(__name__)
 
 
-def _fetch_ic(x0, dx0, system, func, eps=0.001):
+def _fetch_ic(x0, dx0, system, func, t0, eps=0.001):
     if isinstance(x0, list):
         assert len(x0) == len(system.state_vars)
         X0 = np.array(x0, dtype=np.float64)
@@ -22,7 +22,7 @@ def _fetch_ic(x0, dx0, system, func, eps=0.001):
             [np.NaN for _ in system.state_vars], dtype=np.float64
         )
         for k, v in x0.items():
-            *_ , idx = str(k).split('_')
+            _, idx = str(k).split('_')
             idx = int(idx)
             X0[idx] = v
     elif isinstance(x0, (int, float, complex)) and len(system.state_vars) == 1:
@@ -39,7 +39,12 @@ def _fetch_ic(x0, dx0, system, func, eps=0.001):
 
     # if we don't have consistent initial conditions; find them if we can
     # fail if we can't
-    f = lambda y: func(y, X0, 0, 0)
+
+    def f(y):
+        res = np.empty_like(X0)
+        func(t0, X0, y, res)
+        return res
+
     if np.linalg.norm(f(DX0)) > eps:
 
         DX0 = broyden1(f, DX0)
@@ -50,6 +55,7 @@ def _fetch_ic(x0, dx0, system, func, eps=0.001):
 
     return X0, DX0
 
+
 def simulate(system,
              timespan,
              x0,
@@ -101,9 +107,7 @@ def simulate(system,
     Raises:
         ModelException, SolverException
     """
-    from .config import config
-    de = config.de
-    j = config.julia
+
     if system.ports:
         raise ModelException(
             "Cannot Simulate %s: unconnected ports %s",
@@ -112,140 +116,114 @@ def simulate(system,
     if system.control_vars and not control_vars:
         raise ModelException("Control variable not specified")
 
-    tspan = tuple(float(t) for t in timespan)
+    samples = int(1 / dt) + 1
+    t = np.linspace(*timespan, samples)
 
+    res, X = _bondgraph_to_residuals(system, control_vars)
+    X0, DX0 = _fetch_ic(x0, dx0, system, res, t[0])
 
-    func_str, diffs = to_julia_function_string(system, control_vars)
-    func = j.eval(func_str)
-    X0, DX0 = _fetch_ic(x0, dx0, system, func)
+    solver_name = 'ida'
+    dae_solver = dae(solver_name, res)
+    sol = dae_solver.solve(t, X0, DX0)
 
+    return t.reshape((samples, 1)), np.transpose(sol.values.y).T
 
-    problem = de.DAEProblem(func, DX0, X0, tspan, differential_vars=diffs)
 
-    sol = de.solve(problem, dense=True, saveat=float(dt))
+def _to_function(string, X, DX, substitutions):
+    f = sp.sympify(string).subs(substitutions)
 
-    if sol.retcode not in ("Default", "Success"):
-        raise SolverException("Integration error: Solver returned %s "
-                              % sol.retcode, sol)
+    f_n = sp.lambdify((sp.S('t'), X, DX), f, "numpy")
+    return f_n
 
-    t = np.transpose(sol.t)
-
-    return np.resize(t, (len(t), 1)), np.transpose(sol.u).T
-
-
-def to_julia_function_string(model, control_vars=None, in_place=False):
-    """
-    Produces a Julia function string from the given model.
-
-    We expect that that control_vars is a dict with the same keys,
-    or list of the same size, as the model.control_vars
-
-    Args:
-        model:
-        control_vars:
-        in_place:
-
-    Returns:
-        (string, list)
-        A string containing the function definition, and a list of bools
-        identifing which variables contain derivatives.
-    """
 
+def _bondgraph_to_residuals(model, control_vars=None):
     dX = sp.IndexedBase('dX')
     X = sp.IndexedBase('X')
+    U = sp.IndexedBase('U')
     x_subs = []
     dx_subs = []
+    u_subs = []
+    u_func = []
+    n = len(model.state_vars)
+    m = 0
 
     for i, x in enumerate(model.state_vars):
-        x_subs.append((x, X[i+1]))
-        dx_subs.append((sp.S(f'dx_{i}'), dX[i+1]))
-
-
-    cv_strings, dcv_strings = _generate_control_strings(
-        list(model.control_vars.keys()),
-        control_vars,
-        x_subs,
-        dx_subs
-    )
-
-    differential_vars = []
-    subs = x_subs + dx_subs
-
-    function_header = ""
-    function_body = ""
-    function_footer = ""
-
-    if in_place:
-        function_header = "function f(res, dX, X, p, t)\n"
-    else:
-        k = len(model.constitutive_relations)
-        function_header = "function f(dX, X, p, t)\n"
-        function_header += f"    res = zeros({k})\n"
-
-    for cv, dcv in zip(cv_strings, dcv_strings):
-        function_header += cv
-        if dcv:
-            function_header += dcv
-
-    assert model.constitutive_relations
-
-    for relation in model.constitutive_relations:
-
-        eqn_str = str(relation.subs(subs))
-        eqn_str = eqn_str.replace('**', '^')
-        if 'dX' in eqn_str:
-            differential_vars.append(True)
+        x_subs.append((x, X[i]))
+        dx_subs.append((sp.S(f'dx_{i}'), dX[i]))
+
+    if len(model.control_vars) > 0:
+        u_func_dict = {}
+        u_constants = {}
+        if isinstance(control_vars, list):
+            u_func_dict.update({
+                i: f for i, f in enumerate(control_vars)}
+            )
+        elif isinstance(control_vars, dict):
+            u_func_dict.update({
+                int(v[2:]): f for v, f in control_vars.items()
+            })
+        elif len(model.control_vars) == 1:
+            u_func_dict[0] = control_vars
         else:
-            differential_vars.append(False)
-
-        function_body += f"    res[{len(differential_vars)}] = {eqn_str}\n"
-
-    if in_place:
-        function_footer += "end\n"
-    else:
-        function_footer += "    return res\n end\n"
-
-    out_str = function_header + function_body + function_footer
-
-    return (out_str,
-            differential_vars)
+            raise TypeError(f"Control argument {control_vars} is invalid")
+
+        test_x = np.zeros(shape=(n,), dtype=np.float32)
+        for idx, f in u_func_dict.items():
+            try:
+                if isinstance(f, str):
+                    f = _to_function(f, X, dX, dx_subs + x_subs)
+                    u_func_dict[idx] = f
+                if isinstance(f, (float, int, sp.Number)):
+                    u_constants[idx] = f
+                if n == 1:
+                    r = f(0, 0, 0)
+                else:
+                    r = f(0, test_x, test_x)
+                assert isinstance(r, (float, int, sp.Number)
+                                  ), "Invalid output from control"
+            except Exception as ex:
+                message = f"Invalid control function for var: {idx}.\n " \
+                    "Control functions should be of the form:\n" \
+                    f"{idx} = f(t, x, dx/dt)"
+                raise ModelException(message)
+
+        for i, u in enumerate(model.control_vars):
+            if i in u_constants:
+                u_subs.append((u, u_constants[i]))
+                continue
+            u_subs.append((u, U[m]))
+            try:
+                u_func.append(u_func_dict[i])
+            except KeyError:
+                raise ModelException(f"Control variable {u} must be specified")
+            m += 1
+    rels = [r.subs(dx_subs).subs(x_subs).subs(u_subs)
+            for r in model.constitutive_relations]
+
+    if len(rels) != n:
+        raise ModelException(
+            "Model simplification error: system is under-determined")
 
+    Fsym = sp.symarray('F', shape=n)
+    for i, r in enumerate(rels):
+        Fsym[i] = r
 
-def _generate_control_strings(cv, cv_substitutions, x_subs, dx_subs):
+    t = sp.S('t')
+    _r = np.empty(shape=(n,), dtype=np.float64)
+    if not u_func:
+        F = sp.lambdify((t, X, dX), Fsym)
 
-    if isinstance(cv_substitutions, dict):
-        pairs = [(cv_i, cv_substitutions[cv_i]) for cv_i in cv]
-    elif isinstance(cv_substitutions, list):
-        pairs = list(zip(cv, cv_substitutions))
-    elif len(cv) == 1 and (isinstance(cv_substitutions, (str, sp.Expr))):
-        pairs = [cv[0], cv_substitutions]
-    elif not cv and not cv_substitutions:
-        return [], []
+        def residual(_t, _x, _dx, _res):
+            _r = F(_t, _x, _dx)
+            for i in range(n):
+                _res[i] = _r[i]
     else:
-        raise NotImplementedError(
-            f"Could not substitute {cv_substitutions} into {cv}"
-        )
-
-    cv_strings = []
-    dcv_strings = []
-    subs = x_subs + dx_subs
-    partial_pairs = [(X_i,DX_i) for (_,X_i), (_,DX_i) in zip(x_subs, dx_subs)]
-    for var, val in pairs:
-        try:
-
-            u_i = sp.sympify(val).subs(subs)
-            du_i = u_i.diff(sp.S('t')) + sum([
-                u_i.diff(X)*DX for X,DX in partial_pairs
-            ])
-            cv_strings.append(f"    {str(var)} = {u_i}\n")
-            dcv_strings.append(f"    d{str(var)} = {du_i}\n")
-        except SympifyError as ex:
-            s = val
-
-            for x, X in reversed(x_subs):
-                s = s.replace(str(x), str(X))
-
-            cv_strings.append(f"    {str(var)} = {s}\n")
-            dcv_strings.append(None)
-
-    return cv_strings, dcv_strings
+        _u = np.empty(shape=(m,), dtype=np.float64)
+        Fsym_u = sp.lambdify((t, X, dX, U), Fsym)
+
+        def residual(_t, _x, _dx, _res):
+            _u = [u_f(_t, _x, _dx) for u_f in u_func]
+            _r = Fsym_u(_t, _x, _dx, _u)
+            for i in range(n):
+                _res[i] = _r[i]
+    return residual, X
diff --git a/BondGraphTools/version.py b/BondGraphTools/version.py
index 771bc6e..3d26edf 100644
--- a/BondGraphTools/version.py
+++ b/BondGraphTools/version.py
@@ -1,1 +1,1 @@
-__version__ = "0.3.9"
+__version__ = "0.4.1"
diff --git a/BondGraphTools/view.py b/BondGraphTools/view.py
index 94ad4ca..b8f4793 100644
--- a/BondGraphTools/view.py
+++ b/BondGraphTools/view.py
@@ -48,7 +48,7 @@ def _build_graph(system):
 
     try:
         comp_map = {comp: i for i, comp in enumerate(system.components)}
-        graph = dok_matrix((len(comp_map), len(comp_map)), dtype=np.int)
+        graph = dok_matrix((len(comp_map), len(comp_map)), dtype=int)
         for (c1, _), (c2, _) in system.bonds:
             graph[(comp_map[c1], comp_map[c2])] = 1
             graph[(comp_map[c2], comp_map[c1])] = 1
@@ -65,7 +65,7 @@ def _networkx_layout(graph):
     nx_graph = nx.Graph(graph)
     #layout = nx.spring_layout(nx_graph, k=1)
     layout = nx.kamada_kawai_layout(nx_graph, scale=20)
-    pos = [(pair[0],pair[1]) for pair in list(layout.values())]
+    pos = [(pair[0], pair[1]) for pair in list(layout.values())]
 
     return pos
 
@@ -84,11 +84,11 @@ def __init__(self, ax, string, pos, dir, text_dict):
 
         self.x, self.y = pos
         if dir == 'top':
-            self.y += self.height/2
+            self.y += self.height / 2
         elif dir == 'bottom':
             self.y -= self.height / 2
         elif dir == 'right':
-            self.x += self.width/2
+            self.x += self.width / 2
         else:
             self.x -= self.width / 2
 
@@ -108,10 +108,10 @@ def __init__(self, node):
         self.height = 0.1
         self._text = None
         self.ports = {
-            'top':[],
+            'top': [],
             'right': [],
             'bottom': [],
-            'left':[]
+            'left': []
         }
 
     @property
@@ -143,36 +143,36 @@ def axes(self, ax):
 
     def add_port(self, string, dir):
 
-        dx,dy = dir
+        dx, dy = dir
         text_dict = {
             'size': FONT_SM
         }
 
         if dy > abs(dx):
             text_dict.update({
-                'xytext':(self.x, self.y+self.height/2),
-                'horizontalalignment':'center',
+                'xytext': (self.x, self.y + self.height / 2),
+                'horizontalalignment': 'center',
                 'verticalalignment': 'bottom'
-                }
+            }
             )
             dir = 'top'
         elif -dy > abs(dx):
 
             text_dict.update({
-                'xytext':(self.x, self.y-self.height/2),
-                'horizontalalignment':'center',
+                'xytext': (self.x, self.y - self.height / 2),
+                'horizontalalignment': 'center',
                 'verticalalignment': 'top'
-                }
+            }
             )
-            dir= 'bottom'
+            dir = 'bottom'
 
         elif dx >= abs(dy):
 
             text_dict.update({
-                'xytext':(self.x+self.width/2, self.y),
+                'xytext': (self.x + self.width / 2, self.y),
                 'horizontalalignment': 'left',
                 'verticalalignment': 'center'
-                }
+            }
             )
             dir = 'right'
         else:
@@ -184,7 +184,7 @@ def add_port(self, string, dir):
             )
             dir = 'left'
 
-        port = PortGlyph(self.axes, string, self.pos,dir, text_dict)
+        port = PortGlyph(self.axes, string, self.pos, dir, text_dict)
         self.ports[dir] = port
 
         return port
@@ -198,7 +198,7 @@ class BondView(Line2D):
     def __init__(self, port_1, port_2, *args, **kwargs):
         self.port_1 = port_1
         self.port_2 = port_2
-        super().__init__([],[],*args, **kwargs)
+        super().__init__([], [], *args, **kwargs)
 
     def calc_lines(self):
         x1, y1 = self.port_1.pos
@@ -209,12 +209,12 @@ def calc_lines(self):
 
         dx = x2 - x1
         dy = y2 - y1
-        x1 += r1*dx
-        y1 += r1*dy
-        x2 -= r2*dx
+        x1 += r1 * dx
+        y1 += r1 * dy
+        x2 -= r2 * dx
         y2 -= r2 * dy
 
-        lx, ly = x2-x1, y2-y1
+        lx, ly = x2 - x1, y2 - y1
         L = np.sqrt((lx)**2 + (ly)**2)
 
         if not self.shortest_bond or self.shortest_bond > L:
@@ -223,12 +223,12 @@ def calc_lines(self):
         lx /= L
         ly /= L
 
-        headlength = self.shortest_bond/5
+        headlength = self.shortest_bond / 5
 
         vect = np.array((lx, ly))
 
         assert abs(np.linalg.norm(vect) - 1) < 0.01
-        x3, y3 = headlength *self.R.dot(vect) + (x2, y2)
+        x3, y3 = headlength * self.R.dot(vect) + (x2, y2)
 
         self.set_xdata([x1, x2, x3])
         self.set_ydata([y1, y2, y3])
@@ -248,26 +248,26 @@ def draw(self, ax, layout=_networkx_layout):
         ax.get_yaxis().set_visible(False)
         ax.get_xaxis().set_visible(False)
 
-        for component, (x,y) in zip(self._node.components,
-                                                  points):
+        for component, (x, y) in zip(self._node.components,
+                                     points):
             x_min = min(x, x_min)
             x_max = max(x, x_max)
             y_min = min(y, y_min)
             y_max = max(y, y_max)
 
-            component.view.pos = (x,y)
-            if component.metamodel not in {'0','1'}:
+            component.view.pos = (x, y)
+            if component.metamodel not in {'0', '1'}:
                 try:
-                    component.view.string = "\mathbf{{{t}}}: {n}".format(
+                    component.view.string = r"\mathbf{{{t}}}: {n}".format(
                         t=component.metamodel, n=component.name)
-                except:
+                except BaseException:
                     component.view.string = "{t}: {n}".format(
                         t=component.metamodel, n=component.name)
             else:
                 try:
-                    component.view.string = "\mathbf{{{t}}}".format(
+                    component.view.string = r"\mathbf{{{t}}}".format(
                         t=component.metamodel)
-                except:
+                except BaseException:
                     component.view.string = "{t}".format(
                         t=component.metamodel, n=component.name)
             component.view.axes = ax
@@ -302,25 +302,25 @@ def draw(self, ax, layout=_networkx_layout):
         width = abs(x_max - x_min)
         height = abs(y_min - y_max)
         tweak = 0.1
-        ax.axis([x_min - tweak*width,
-                 x_max + tweak*width,
-                 y_min - tweak*height,
-                 y_max + tweak*height])
+        ax.axis([x_min - tweak * width,
+                 x_max + tweak * width,
+                 y_min - tweak * height,
+                 y_max + tweak * height])
 
 
 def find_renderer(fig):
 
     if hasattr(fig.canvas, "get_renderer"):
-        #Some backends, such as TkAgg, have the get_renderer method, which
-        #makes this easy.
+        # Some backends, such as TkAgg, have the get_renderer method, which
+        # makes this easy.
         renderer = fig.canvas.get_renderer()
     else:
-        #Other backends do not have the get_renderer method, so we have a work
-        #around to find the renderer.  Print the figure to a temporary file
-        #object, and then grab the renderer that was used.
-        #(I stole this trick from the matplotlib backend_bases.py
-        #print_figure() method.)
+        # Other backends do not have the get_renderer method, so we have a work
+        # around to find the renderer.  Print the figure to a temporary file
+        # object, and then grab the renderer that was used.
+        # (I stole this trick from the matplotlib backend_bases.py
+        # print_figure() method.)
         import io
         fig.canvas.print_pdf(io.BytesIO())
         renderer = fig._cachedRenderer
-    return(renderer)
\ No newline at end of file
+    return(renderer)
diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 0000000..488a969
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,15 @@
+FROM ubuntu:focal
+
+ENV TZ=Australia/Melbourne
+RUN ln -snf /usr/share/zoneinfo/$TZ /etc/localtime && echo $TZ > /etc/timezone
+
+RUN apt-get update && apt-get install -y gfortran libblas-dev liblapack-dev git build-essential cmake
+
+RUN git clone https://github.com/LLNL/sundials.git; \
+    mkdir /sundials/build;\
+    cd /sundials/build; \
+    cmake -DLAPACK_ENABLE=ON ..; \
+    make && make install
+
+RUN apt-get install -y python3-dev python3-pip
+RUN python3 -m pip install bondgraphtools
diff --git a/README.md b/README.md
index cd83590..d8a283d 100644
--- a/README.md
+++ b/README.md
@@ -15,26 +15,23 @@ https://bondgraphtools.readthedocs.io/
 
 ### Dependencies
 
-BondGraph requires:
-- python 3.6
-- julia 0.6.4
+BondGraphTools requires:
+- python 3.7
+- LaTeX 
+- gFortran, BLAS, LAPACK and SUNDIALS (for the numerical solver)
 
 Python dependencies:
-- sympy>=1.1.1
-- numpy>=1.14
-- scipy>=1.0.1
-- matplotlib>=2.2.2
-- julia>=0.1.5
-- diffeqpy>=0.4
-
-Julia dependencies:
- - PyCall
- - DifferentialEquations.jl
-
-### Instructions:
-1. Install python > 3.6 for your operating system.
-2. Install Julia 0.6.4 (https://julialang.org/downloads/) for your operating
- system. _Julia 0.7 and 1.0 are not yet supported_
-3. Make sure Julia 0.6.4 is in your os path. (test this by running `julia -v`)
-4. Install using PyPI; `pip install BondGraphTools`
-
+- sympy
+- numpy
+- scipy
+- matplotlib
+- scikits.odes 
+
+## Instructions:
+
+#### Recommended installation using Conda
+1. Install [Anaconda](https://anaconda.org/) (or miniconda).
+2. Install dependencies from conda-forge: liblas, liblapack, sundials using
+   `conda install -c conda-forge libblas liblapack sundials fortran-compiler`
+3. Install LaTeX via `conda install -c conda-forge miktex` (windows) or `conda install -c conda-forge texlive-core` (osx/linux)
+4. Install the package using `pip install bondgraphtools`
diff --git a/docs/api.rst b/docs/api.rst
index 8482b4a..e8826b5 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -51,7 +51,7 @@ API Reference
 :mod:`BondGraphTools.sim_tools`
 -------------------------------
 .. automodule:: BondGraphTools.sim_tools
-    :members: to_julia_function_string
+    :members: simulate
     :show-inheritance:
 
 :mod:`BondGraphTools.reaction_builder`
diff --git a/docs/guides.rst b/docs/guides.rst
deleted file mode 100644
index ad70ca3..0000000
--- a/docs/guides.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-How-To Guides
-=============
-
-.. toctree::
-   :maxdepth: 1
-   :glob:
-
-   guides/*
diff --git a/docs/guides/installing-on-linux.rst b/docs/guides/installing-on-linux.rst
deleted file mode 100644
index 2b93b3b..0000000
--- a/docs/guides/installing-on-linux.rst
+++ /dev/null
@@ -1,93 +0,0 @@
-Installing BondGraphTools on Linux
-===================================
-
-A working installation of both Python and Julia is required to use all of
-BondGraphTools features.
-In this guide we walk step-by-step through the recommended method of
-installation.
-
-Commands in this section beginning with `$` are entered at the terminal prompt.
-Commands beginning with `>>>` are enter at the Python interpreter prompt.
-
-We assume Ubuntu 16.04 or greater. BondGraphTools _should_ work on all
-distributions but it is only tested and supported on Ubuntu 16.04 and up.
-
-Step 1: Install Python 3.6 or greater
--------------------------------------
-
-Ubuntu 17.10 or Greater
-+++++++++++++++++++++++
-Skip this step as Python 3.6 is already installed.
-
-Ubuntu 14.04 Only
-+++++++++++++++++
-A newer version of the C++ standard library is required. First, download it a
-working directory with::
-
-  $ wget -q -O libstdc++6 http://security.ubuntu.com/ubuntu/pool/main/g/gcc-5/libstdc++6_5.4.0-6ubuntu1~16.04.10_amd64.deb
-
-then install using::
-
-    $ sudo dpkg --force-all -i libstdc++6
-
-Ubuntu 14.04 and 16.04
-++++++++++++++++++++++
-Add the 'deadsnakes' repository to the Ubuntu package manager::
-
-    $ sudo add-apt-repository ppa:deadsnakes/ppa
-
-Ubuntu 14.04 to 16.10
-+++++++++++++++++++++
-Update the package indices and install Python 3.6::
-
-    $ sudo apt update
-    $ sudo apt install python36 python36-pip
-
-Step 2: Install Julia
----------------------
-We must now install Julia 0.6.
-Our examples follow the Julia platform instructions https://julialang.org/downloads/platform.html#generic-binaries .
-
-1. Download Julia 0.6.4 generic binaries from https://julialang.org/downloads/oldreleases.html .
-   For example using curl to save the archive to your downloads::
-
-    $ curl -L https://julialang-s3.julialang.org/bin/linux/x64/0.6/julia-0.6.4-linux-x86_64.tar.gz > ~/Downloads/julia-0.6.4-linux-x86_64.tar.gz
-2. Extract the archive `.tar.gz` to the desired directory. For example, here
-   we first create the target directory, then extract the tarball::
-
-    $ sudo mkdir /usr/local/opt/julia
-    $ sudo tar xf ~/Downloads/julia-0.6.4-linux-x86_64.tar.gz --directory /usr/local/opt/julia --strip-components=1
-
-3. Add julia to the PATH, for example by adding a symbolic link to the local bin
-   folder::
-
-    $ sudo ln -s /usr/local/opt/julia/bin/julia /usr/local/bin/julia
-
-4. Verify that this has worked by typing `julia -v', which should result in
-   `julia version 0.6.4`.
-
-Step 3: Install BondGraphTools
-------------------------------
-Use python to install BondGraphTools with::
-
-    $ pip3 install BondGraphTools
-
-or alternatively::
-
-    $ python3.6 -m pip install BondGraphTools
-
-It is optional but highly recommended to install the Jupyter packages via::
-
-    $ pip3 install Jupyter
-
-Step 4: Julia dependencies.
----------------------------
-`BondGraphTools` will automatically download and install Julia dependencies the first time
-a user attempts to run a simulation.
-
-This can be manually triggered via the following commands in the python
-interpreter::
-
-    >>> from bondgraphtools.config import config
-    >>> config.install_dependencies()
-
diff --git a/docs/guides/installing-on-osx.rst b/docs/guides/installing-on-osx.rst
deleted file mode 100644
index 2816364..0000000
--- a/docs/guides/installing-on-osx.rst
+++ /dev/null
@@ -1,60 +0,0 @@
-Installing BondGraphTools on Mac OSX
-====================================
-
-A working installation of both Python and Julia is required to use all of BondGraphTools features.
-In this guide we walk step-by-step through the recommended method of installation.
-
-Step 1: Instal Homebrew and Python
-----------------------------------
-First we must install python 3. We recommend doing this via Homebrew_.
-
-1. Open a terminal instance.
-2. Paste the following command into the terminal instance::
-
-    $ /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
-
-3. Once Homebrew has finished, install python by typing the following command into the terminal::
-
-    $ brew install python
-
-Python 3.7 is now installed.
-
-Step 2: Install Julia
----------------------
-We must now install Julia 0.6.
-
-1. Download Julia 0.6.4 from https://julialang.org/downloads/oldreleases.html
-2. Once the download is complete open the file. It will contain a link to th Applications folder, and the Julia-0.6
-   package. Drag the julia package into the applications folder. Once this is done, julia has been installed and you
-   can remove the downloaded file.
-3. Add julia to your profile path by typing the following command into the terminal::
-
-    $ echo 'export PATH=/Applications/Julia-0.6.app/Contents/Resources/julia/bin/:$PATH' >> ~/.bash_profile
-
-Step 3: Install BondGraphTools
-------------------------------
-
-1. Install BondGraphTools from PyPI by typing the following command into the terminal::
-
-    $ pip3 install BondGraphTools
-
-2. Install the recommended additional packages (jupyter) by typing the following command into the terminal::
-
-    $ pip3 install jupyter
-
-`BondGraphTools` is now installed.
-
-Step 4: Julia dependencies.
----------------------------
-`BondGraphTools` will automatically download and install Julia dependencies the first time
-a user attempts to run a simulation.
-
-This can be manually triggered via the following command::
-
-    >>> from bondgraphtools.config import config
-    >>> config.install_dependencies()
-
-.. _Homebrew: https://brew.sh/
-
-
-
diff --git a/docs/guides/troubleshooting.rst b/docs/guides/troubleshooting.rst
deleted file mode 100644
index 1a5743b..0000000
--- a/docs/guides/troubleshooting.rst
+++ /dev/null
@@ -1,18 +0,0 @@
-Troubleshooting
-===============
-
-Julia Errors
-------------
-If you encounter julia errors, you may need to rebuild pycall, or re-install your julia libraries.
-
-1. To rebuild pycall and related libraries enter the following commands in your python interpreter::
-
-        from bondgraphtools.config import config
-        config.install_dependencies(rebuilt=True)
-
-2. If that fails, you may need to reset your julia library. **Warning** this will delete any
-   additional libraries you may have installed. This can be done by removing the directory '~/.julia/lib'
-   in Linux/MacOSX, or 'C:\\users\\<username>\\.julia\\lib' in windows.
-
-
-
diff --git a/docs/index.rst b/docs/index.rst
index d38f7c5..549a42a 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -3,7 +3,6 @@
    :maxdepth: 1
 
    tutorials
-   guides
    discussion
    api
 
@@ -19,9 +18,8 @@ BondGraphTools is intended to be used by:
    * engineers; to quickly build and simulate physical systems
    * mathematicians; to perform model reduction and analysis
 
-BondGraphTools is built upon the scientific python stack; numpy, scipy, sympy
-and matplotlib, and relies on Julia and DifferentialEquations.jl to handle
-numnerical simulation.
+BondGraphTools is built upon the scientific python stack; numpy, scipy, sympy,
+scikit.odes and matplotlib
 
 How to read the docs
 ^^^^^^^^^^^^^^^^^^^^
@@ -38,21 +36,16 @@ Getting Started
 Requirements
 ------------
 
-* Python 3.6 or above.
-* Julia 0.6. *Note: Julia 0.7 and above are not yet supported*
+* Python 3.7 or above.
 
 Installation
 ------------
 
-1. Install Julia 0.6.4 which can be downloaded from here_.
-2. Make sure Julia 0.6.4 is in your path variables. You can test this by typing
-   ``julia -v`` at the command prompt. If this returns ``julia version 0.6.4``
-   then julia is in your path variables.
+1. Install Python 3.7 or greater.
+2. Install Sundials and ``scikits.odes`` [dependencies](https://scikits-odes.readthedocs.io/en/stable/installation.html)
 3. Install BondGraphTools via pip: ``pip install BondGraphTools``
 
-BondGraphTools is now installed and will download and install the required
-julia packages the first time you import it.
-
+BondGraphTools is now installed.
 
 Usage
 -----
@@ -62,9 +55,6 @@ script by using::
 
     import BondGraphTools
 
-Users need an internet connection during the first load as the it will download
-the required packages from the julia package repository.
-
 Api reference for the package and it's contents can be accessed using the help
 command::
 
@@ -78,6 +68,5 @@ If you wish to help out, please visit the project github_.
 Bibliography
 ^^^^^^^^^^^^
 .. _github: https://github.com/BondGraphTools/BondGraphTools
-.. _here: https://julialang.org/downloads/oldreleases.html
 .. [Gaw2007] https://ieeexplore.ieee.org/document/4140745
 .. [Gaw1996] http://www.gawthrop.net/Books/GawSmi96.pdf
diff --git a/docs/tutorials/RC.rst b/docs/tutorials/RC.rst
index e996a59..5718351 100644
--- a/docs/tutorials/RC.rst
+++ b/docs/tutorials/RC.rst
@@ -6,14 +6,14 @@ Tutorial: Driven Filter Circuit
 +----------------+------------------------------------------------------------+
 | Difficulty:    | Beginner.                                                  |
 +----------------+------------------------------------------------------------+
-| Requirement:   | `BondGraphTools`, `jupyter`.                               |
+| Requirement:   | ``BondGraphTools``, Jupyter.                               |
 +----------------+------------------------------------------------------------+
-| How to follow: | Enter each block of code in consecutive cells in a jupyter |
+| How to follow: | Enter each block of code in consecutive cells in a Jupyter |
 |                | notebook.                                                  |
 +----------------+------------------------------------------------------------+
 
 In part 1 of this tutorial we will demonstrate how to build and connect models
-using `BondGraphTools` by constructing a simple passive filter.
+using ``BondGraphTools`` by constructing a simple passive filter.
 Part 2 introduces control sources, and provides examples of how one can perform
 parameter sweeps or input comparisons.
 
@@ -21,7 +21,7 @@ parameter sweeps or input comparisons.
 Part 1: Basic Use
 -----------------
 
-First, import `BondGraphTools` and create a new model with the name "RC"::
+First, import ``BondGraphTools`` and create a new model with the name "RC"::
 
     import BondGraphTools as bgt
     model = bgt.new(name="RC")
@@ -58,14 +58,14 @@ which produces a sketch of the network topology.
 
 
 To demonstrate that the isolated system is behaving correctly, we simulate
-from the initial where the C component has $x_0=1$ and run the simulation over
-the  time interval $(0,5)$. This results in a vector $t$ of time of time points
-and a corresponding vector $x$ of data points which can then be plotted against
-eachother with `matplotlib` ::
+from the initial where the C component has :math:`x_0=1` and run the simulation over
+the  time interval :math:`(0,5)`. This results in a vector :math:`t` of time points
+and a corresponding vector :math:`x` of data points which can then be plotted against
+each other with ``matplotlib``::
 
     timespan = [0, 5]
     x0 = [1]
-    t, x = simulate(model, timespan=timespan, x0=x0)
+    t, x = bgt.simulate(model, timespan=timespan, x0=x0)
     from matplotlib.pyplot import plot
     fig = plot(t,x)
 
@@ -77,47 +77,47 @@ Part 2: Control
 ---------------
 We wish to see how this filter responds to input.
 Add flow source by creating a new Sf component, adding to the model, and
-connecting it to the common voltage law.::
+connecting it to the common voltage law::
 
     Sf = bgt.new('Sf')
     bgt.add(model, Sf)
     bgt.connect(Sf, zero_law)
 
-The model should now look something like this ::
+The model should now look something like this::
 
-    draw(model)
+    bgt.draw(model)
 
 .. figure:: images/RC_3.svg
     :scale: 50 %
     :align: center
 
-The model also now has associated with it a control variable `u_0`.
-Control variables can be listed via the attribute `model.control_vars` and we
+The model also now has associated with it a control variable ``u_0``.
+Control variables can be listed via the attribute ``model.control_vars`` and we
 can observe the constitutive relations, which give the implicit equations of
 motion for the system in sympy form::
 
     model.constitutive_relations
     # returns [dx_0 - u_0 + x_0]
 
-where `x_0` and `dx_0` are the state variable and it's derivative. One can
+where ``x_0`` and ``dx_0`` are the state variable and it's derivative. One can
 identify where that state variable came from via::
 
     model.state_vars
     # returns {'x_0': (C: C1, 'q_0)}
 
-Here `C: C1` is a reference to the `C` object itself.
+Here ``C: C1`` is a reference to the ``C`` object itself.
 
 
 Part 3: Simulations
 -------------------
 We will now run various simulations.
 
-Firstly, we simulate with constant effort by passing the control law $u_0=2$
+Firstly, we simulate with constant effort by passing the control law :math:`u_0=2`
 to the solver and plotting the results::
 
     timespan = [0, 5]
     x0 = [1]
-    t, x = simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':2})
+    t, x = bgt.simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':2})
     plot(t,x)
 
 .. figure:: images/RC_4.svg
@@ -125,21 +125,20 @@ to the solver and plotting the results::
     :align: center
 
 Time dependent control laws can be specified as string. In this case we
-consider the response to a $\pi^{-1}$Hz sine wave.::
+consider the response to a :math:`\pi^{-1}` Hz sine wave::
 
-    t, x = simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':'sin(2*t)'})
+    t, x = bgt.simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':'sin(2*t)'})
     plot(t,x)
 
-
 .. figure:: images/RC_5.svg
     :scale: 50 %
     :align: center
 
 One can also consider the impulse response of by applying a step function
-input ot the control law.::
+input ot the control law::
 
     step_fn = 't < 1 ? 1 : 0' # if t < 0 then 1 else 0
-    t, x = simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':step_fn})
+    t, x = bgt.simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':step_fn})
     plot(t,x)
 
 .. figure:: images/RC_6.svg
@@ -147,12 +146,12 @@ input ot the control law.::
     :align: center
 
 Finally we run a sequence of simulations where a new control law is generated
-based on the loop iteration.
+based on the loop iteration::
 
     fig = plt.figure()
     for i in range(4):
         func_text = "cos({i}t)".format(i=i)
-        t_i, x_i = simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':func_text})
+        t_i, x_i = bgt.simulate(model, timespan=timespan, x0=x0, control_vars={'u_0':func_text})
         plot(t_i,x_i)
 
 .. figure:: images/RC_7.svg
diff --git a/requirements.txt b/requirements.txt
index 5747e06..7db2644 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,10 +1,9 @@
-sympy>=1.1
-numpy==1.14.5
-scipy>=1.0.1
-matplotlib==2.2.2
-networkx>=2.1
-requests>=2.18.4
-julia==0.1.5
-diffeqpy==0.4
-pyYAML>=3.0
-ordered-set>=3.0.0
\ No newline at end of file
+numpy
+sympy
+scipy
+matplotlib
+networkx
+requests
+pyYAML
+ordered-set
+scikits.odes
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 9434eb8..220611d 100644
--- a/setup.py
+++ b/setup.py
@@ -24,6 +24,8 @@
         'Development Status :: 3 - Alpha',
         "Programming Language :: Python :: 3.6",
         "Programming Language :: Python :: 3.7",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
         "Topic :: Scientific/Engineering :: Mathematics",
         "Topic :: Scientific/Engineering",
         "License :: OSI Approved :: Apache Software License",


Test patch:
diff --git a/test/atpase.py b/test/atpase.py
index f08ba58..03db749 100644
--- a/test/atpase.py
+++ b/test/atpase.py
@@ -1,11 +1,13 @@
+import BondGraphTools as bgt
 import pathlib
 import sys
 import os
 path = os.path.relpath(path=".")
 sys.path.insert(0, os.path.abspath(path))
 
-import BondGraphTools as bgt
 log_file = "profile.log"
+
+
 def atpase():
 
     # Reaction Conductance in fmol/s
@@ -93,10 +95,10 @@ def atpase():
     model = bgt.new(name="ATPase")
     # add all the pumps
     for pump, affinity in enumerate(pump_affinities):
-        parameters = {'k':affinity,
-                      'R':Gas_Constant,
-                      'T':Temperature}
-        name = "P{}".format(pump+1)
+        parameters = {'k': affinity,
+                      'R': Gas_Constant,
+                      'T': Temperature}
+        name = "P{}".format(pump + 1)
         Ce = bgt.new(component="Ce",
                      name=name,
                      library="BioChem",
@@ -108,9 +110,9 @@ def atpase():
 
     # add all the species
     for species, affinity in species_affinities.items():
-        parameters = {'k':affinity,
-                      'R':Gas_Constant,
-                      'T':Temperature}
+        parameters = {'k': affinity,
+                      'R': Gas_Constant,
+                      'T': Temperature}
         Ce = bgt.new(component="Ce",
                      name=species,
                      library="BioChem",
@@ -120,17 +122,19 @@ def atpase():
         model.add(Ce, junction)
         model.connect(Ce, junction)
 
-    #add the electrical components
-    c_mem = bgt.new('C',name="mem",value={'C':membrane_capacitance})
+    # add the electrical components
+    c_mem = bgt.new('C', name="mem", value={'C': membrane_capacitance})
     c_junction = bgt.new("0", name="mem")
     model.add(c_mem, c_junction)
     model.connect(c_mem, c_junction)
 
-    zF_5 = bgt.new('TF',name='zF_5', value={'r':r5_charge})
-    zF_8 = bgt.new('TF',name='zF_8', value={'r':r8_charge})
+    zF_5 = bgt.new('TF', name='zF_5', value={'r': r5_charge})
+    zF_8 = bgt.new('TF', name='zF_8', value={'r': r8_charge})
 
     model.add(zF_5, zF_8)
-    port= 1 #  port 0 is the primary winding of a transformer and 1 is the secondary winding.
+    # port 0 is the primary winding of a transformer and 1 is the secondary
+    # winding.
+    port = 1
     zF_5_secondary = (zF_5, port)
     zF_8_secondary = (zF_8, port)
     model.connect(zF_5_secondary, c_junction)
@@ -162,7 +166,7 @@ def atpase():
             for product in products:
                 if product == 'zF_5':
                     junction = (
-                    zF_5, 0)  # the primary winding of the charge transporter
+                        zF_5, 0)  # the primary winding of the charge transporter
                 elif product == 'zF_8':
                     junction = (zF_8, 0)
                 else:
@@ -176,11 +180,11 @@ def atpase():
     from math import log
     # Add Chemostats
     for species, affinity in species_affinities.items():
-        effort = log(affinity*initial_conditions[species])
+        effort = log(affinity * initial_conditions[species])
         Se = bgt.new('Se', value=effort, name=species)
         model.add(Se)
         junction = model.find(name=species, component='0')
         model.connect(Se, junction)
 
     # generate equations
-    model_equations = model.constitutive_relations
\ No newline at end of file
+    model_equations = model.constitutive_relations
diff --git a/test/conftest.py b/test/conftest.py
index 6f9aa51..4af94f0 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -7,6 +7,10 @@ def pytest_addoption(parser):
                      default=False, help="run slow tests")
 
 
+def pytest_configure(config):
+    config.addinivalue_line("markers", "slow: mark tests as slow")
+
+
 def pytest_collection_modifyitems(config, items):
     if config.getoption("--runslow"):
         # --runslow given in cli: do not skip slow tests
@@ -32,4 +36,3 @@ def rlc():
     bgt.connect(c, kvl)
 
     return rlc
-
diff --git a/test/profile_atpase.py b/test/profile_atpase.py
index 6c3bb64..6843774 100644
--- a/test/profile_atpase.py
+++ b/test/profile_atpase.py
@@ -1,3 +1,4 @@
+from pyprof2calltree import convert, visualize
 import logging
 from cProfile import Profile
 from atpase import atpase, bgt
@@ -7,7 +8,7 @@
 profiler = Profile()
 logger = bgt.logger
 logger.setLevel(logging.INFO)
-handler = logging.FileHandler(filename=LOG,mode='w')
+handler = logging.FileHandler(filename=LOG, mode='w')
 logger.addHandler(handler)
 
 try:
@@ -15,8 +16,6 @@
 except KeyboardInterrupt as ex:
     handler.flush()
     raise ex
-from pyprof2calltree import convert, visualize
 
 convert(profiler.getstats(), FILENAME)
 visualize(profiler.getstats())
-
diff --git a/test/test_RLC.py b/test/test_RLC.py
index cb8da54..68deb7b 100644
--- a/test/test_RLC.py
+++ b/test/test_RLC.py
@@ -4,11 +4,11 @@
 import BondGraphTools.sim_tools as sim
 
 
-@pytest.mark.use_fixture("rlc")
 def test_build(rlc):
     assert len(rlc.state_vars) == 2
     assert len(rlc.ports) == 0
 
+
 def test_build_rlc():
     r = bgt.new("R", value=1)
     l = bgt.new("I", value=1)
@@ -23,7 +23,6 @@ def test_build_rlc():
     assert len(kvl.ports) == 3
 
 
-@pytest.mark.use_fixture("rlc")
 def test_build_and_drive(rlc):
     se = bgt.new("Se")
     assert len(se.control_vars) == 1
@@ -46,7 +45,7 @@ def test_symbolic_params():
     c = bgt.new("C", value=sympy.symbols('c'))
     kvl = bgt.new("0", name="kvl")
     rlc = bgt.new()
-    rlc.add([r,l, c , kvl])
+    rlc.add([r, l, c, kvl])
 
     bgt.connect(r, kvl)
     bgt.connect(l, kvl)
@@ -59,14 +58,12 @@ def test_symbolic_params():
     assert params & set(sympy.symbols('r, l, c'))
 
 
-@pytest.mark.use_fixture("rlc")
 def test_rlc_con_rel(rlc):
 
     rel = rlc.constitutive_relations
 
     _, v = rlc.state_vars['x_0']
 
-
     if str(v) != 'q_0':
         eq1 = sympy.sympify("dx_0 - x_1")
         eq2 = sympy.sympify("dx_1 + x_0 + x_1")
@@ -80,21 +77,20 @@ def test_rlc_con_rel(rlc):
     assert "x_0" in rlc.state_vars
     assert "x_1" in rlc.state_vars
 
+
 def test_tf():
     l = bgt.new("I", value=1)
     c = bgt.new("C", value=1)
     tf = bgt.new("TF", value=0.5)
     tflc = bgt.new()
-    tflc.add([tf,l, c])
-
-
+    tflc.add([tf, l, c])
 
     bgt.connect(l, (tf, 1))
     bgt.connect(c, (tf, 0))
 
-    c,m,lp,nlp,const = tflc.system_model()
-    assert nlp.is_zero
-    assert const ==[]
+    c, m, lp, nlp, const = tflc.system_model()
+    assert nlp.is_zero_matrix
+    assert const == []
 
 
 def test_se():
@@ -106,10 +102,10 @@ def test_se():
     assert Se.constitutive_relations == [sympy.sympify("e_0 - 1")]
     bgt.connect(Se, c)
 
-
     assert vc.constitutive_relations == [sympy.sympify("dx_0"),
                                          sympy.sympify("x_0 - 1")]
 
+
 def test_one():
     loop_law = bgt.new('1')
     Se = bgt.new('Se', value=1)
@@ -125,5 +121,3 @@ def test_one():
     assert vc.constitutive_relations == [
         sympy.sympify("dx_0 + x_0 - 1")
     ]
-
-
diff --git a/test/test_algebra.py b/test/test_algebra.py
index 81c01fb..c6b05b8 100644
--- a/test/test_algebra.py
+++ b/test/test_algebra.py
@@ -5,9 +5,10 @@
 import BondGraphTools as bgt
 from BondGraphTools import connect, new, expose
 from BondGraphTools.algebra import extract_coefficients, smith_normal_form, \
-    adjacency_to_dict, augmented_rref,_generate_substitutions,\
+    adjacency_to_dict, augmented_rref, _generate_substitutions,\
     inverse_coord_maps, _generate_cv_substitutions, get_relations_iterator
 
+
 def test_extract_coeffs_lin():
     eqn = sympy.sympify("y -2*x -3")
     local_map = {
@@ -38,8 +39,8 @@ def test_extract_coeffs_nlin():
 
 def test_smith_normal_form():
 
-    m = sympy.SparseMatrix(2,3,{(0,2):2, (1,1):1})
-    mp,_,_ = smith_normal_form(m)
+    m = sympy.SparseMatrix(2, 3, {(0, 2): 2, (1, 1): 1})
+    mp, _, _ = smith_normal_form(m)
     assert mp.shape == (3, 3)
     assert mp[2, 2] != 0
 
@@ -48,19 +49,20 @@ def test_smith_normal_form_2():
     matrix = sympy.eye(3)
     matrix.row_del(1)
 
-    m,_,_ = smith_normal_form(matrix)
+    m, _, _ = smith_normal_form(matrix)
 
-    diff = sympy.Matrix([[0,0,0],[0,1,0], [0,0,0]])
+    diff = sympy.Matrix([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
 
     assert (sympy.eye(3) - m) == diff
 
 
 def test_smith_normal_form_3():
-    m = sympy.SparseMatrix(5,3,{(0,1):1, (1,0):1,
-                                (4,2):1})
-    mp,_,_ = smith_normal_form(m)
+    m = sympy.SparseMatrix(5, 3, {(0, 1): 1, (1, 0): 1,
+                                  (4, 2): 1})
+    mp, _, _ = smith_normal_form(m)
     assert mp.shape == (3, 3)
-    assert (mp - sympy.eye(3)).is_zero
+
+    assert (mp - sympy.eye(3)).is_zero_matrix
 
 
 def test_aug_rref():
@@ -84,47 +86,47 @@ def test_aug_rref():
     assert m != matrix
     assert a != adj
 
-    assert (m - sympy.eye(4)).is_zero
+    assert (m - sympy.eye(4)).is_zero_matrix
 
-    assert (a * matrix - sympy.eye(4)).is_zero
+    assert (a * matrix - sympy.eye(4)).is_zero_matrix
 
 
 def test_augmented_rref():
     M = sympy.Matrix(
         [[1, 0, 1, 0],
          [1, 0, 1, 0],
-         [0, 0, 1 ,0]])
+         [0, 0, 1, 0]])
 
-    A = sympy.Matrix(3,1, list(sympy.symbols('a,a,c')))
+    A = sympy.Matrix(3, 1, list(sympy.symbols('a,a,c')))
     target_A = sympy.Matrix(3, 1,
-                            [sympy.sympify('a - c'), sympy.symbols('c'),0])
+                            [sympy.sympify('a - c'), sympy.symbols('c'), 0])
     Mr, _ = M.rref()
 
-    assert Mr == sympy.Matrix([[1,0,0,0],
-                               [0,0,1,0],
-                               [0,0,0,0]])
+    assert Mr == sympy.Matrix([[1, 0, 0, 0],
+                               [0, 0, 1, 0],
+                               [0, 0, 0, 0]])
 
     MA = M.row_join(A)
     MA = augmented_rref(MA, 1)
-    assert MA[:,:-1]== Mr
-    assert target_A == MA[:,-1:]
+    assert MA[:, :-1] == Mr
+    assert target_A == MA[:, -1:]
 
 
 def test_augmented_rref_2():
     M = sympy.Matrix(
         [[1, 0, 1, 0],
          [1, 0, 1, 0],
-         [0, 0, 1 ,0]])
+         [0, 0, 1, 0]])
 
-    A = sympy.Matrix(3,1, [sympy.symbols('a'),sympy.symbols('a'), 1])
+    A = sympy.Matrix(3, 1, [sympy.symbols('a'), sympy.symbols('a'), 1])
     target_A = sympy.Matrix(3, 1,
                             [sympy.sympify('a - 1'), 1, 0])
     Mr, _ = M.rref()
 
-    assert Mr == sympy.Matrix([[1,0,0,0],
-                               [0,0,1,0],
-                               [0,0,0,0]])
-    assert M.shape == (3,4)
+    assert Mr == sympy.Matrix([[1, 0, 0, 0],
+                               [0, 0, 1, 0],
+                               [0, 0, 0, 0]])
+    assert M.shape == (3, 4)
     assert A.shape == (3, 1)
     MA = M.row_join(A)
     MA = augmented_rref(MA, 1)
@@ -141,6 +143,7 @@ def test_build_relations():
 
     assert set(eqns) == test_eqn
 
+
 def test_zero_junction_relations():
     r = bgt.new("R", value=sympy.symbols('r'))
     l = bgt.new("I", value=sympy.symbols('l'))
@@ -150,7 +153,6 @@ def test_zero_junction_relations():
     rlc = bgt.new()
     rlc.add([c, l, kvl, r])
 
-
     connect(r, kvl)
     connect(l, kvl)
     connect(c, kvl)
@@ -162,7 +164,6 @@ def test_zero_junction_relations():
     assert sympy.sympify("f_0 + f_1 + f_2") in rels
 
 
-@pytest.mark.usefixture('rlc')
 def test_basis_vectors(rlc):
 
     model_basis_vects = set()
@@ -182,8 +183,8 @@ def test_build_junction_dict():
     bg = bgt.new()
     bg.add([c, kvl])
     connect(kvl, c)
-    cp,kp = list(c.ports) + list(kvl.ports)
-    index_map = {cp:0, kp:1}
+    cp, kp = list(c.ports) + list(kvl.ports)
+    index_map = {cp: 0, kp: 1}
     M = adjacency_to_dict(index_map, bg.bonds, offset=1)
     assert M[(0, 1)] == 1
     assert M[(0, 3)] == -1
@@ -204,7 +205,6 @@ def test_build_model_fixed_cap():
     assert test_eqn2 in eqns
 
 
-@pytest.mark.usefixture("rlc")
 def test_rlc_basis_vectors(rlc):
 
     tangent, ports, cv = rlc.basis_vectors
@@ -221,10 +221,10 @@ def test_relations_iter():
     coords = list(sympy.sympify("dx_0,e_0,f_0,x_0, 1"))
     relations = get_relations_iterator(c, mappings, coords)
 
-    d1 = {0:1, 2:-1} #dx/dt - f = 0
-    d2 = {1:1, 3:-1} #e - x = 0$
-    d1m = {0:-1, 2:1} #dx/dt - f = 0
-    d2m = {1:-1, 3:1} #e - x = 0$
+    d1 = {0: 1, 2: -1}  # dx/dt - f = 0
+    d2 = {1: 1, 3: -1}  # e - x = 0$
+    d1m = {0: -1, 2: 1}  # dx/dt - f = 0
+    d2m = {1: -1, 3: 1}  # e - x = 0$
 
     for lin, nlin in relations:
         assert not nlin
@@ -234,7 +234,7 @@ def test_relations_iter():
 def test_relation_iter_twoport():
     tf = bgt.new("TF", value=1)
     tf_port_0, tf_port_1 = list(tf.ports)
-    mappings = ({},{tf_port_0:0, tf_port_1:1},{})
+    mappings = ({}, {tf_port_0: 0, tf_port_1: 1}, {})
     coords = list(sympy.sympify("e_0,f_0,e_1,f_1"))
     # d1 = {0: 1, 2: -1}
     # d2 = {1: 1, 2: 1}
@@ -243,7 +243,6 @@ def test_relation_iter_twoport():
         assert not nlin
 
 
-@pytest.mark.usefixture("rlc")
 def test_interal_basis_vectors(rlc):
     tangent, ports, cv = rlc._build_internal_basis_vectors()
 
@@ -260,7 +259,7 @@ def test_cv_relations():
     bg = bgt.new()
     bg.add([c, se, kcl, r])
 
-    connect(c,(kcl,kcl.non_inverting))
+    connect(c, (kcl, kcl.non_inverting))
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
@@ -286,24 +285,24 @@ def test_parallel_crv_relations():
 def test_generate_subs():
 
     w, x, y, z = sympy.sympify("w,x,y,z")
-    size_tuple  =(0, 2, 0,4 )
+    size_tuple = (0, 2, 0, 4)
     coords = [w, x, y, z]
     #  w + w^2 + x^2
     #  x + 1 + y^2   < should appear in subs as x = -y^2  - 1
     #  y + 1         <                          y = - 1
     #  0 + z^2 + w^2
 
-    L = sympy.SparseMatrix(4, 4, {(0,0): 1,
-                                  (1,1): 1,
-                                  (2,2): 1})
+    L = sympy.SparseMatrix(4, 4, {(0, 0): 1,
+                                  (1, 1): 1,
+                                  (2, 2): 1})
 
     N = sympy.SparseMatrix(4, 1, {(0, 0): w**2 + x**2,
                                   (1, 0): 1 + y**2,
                                   (2, 0): 1})
 
     constraint = [z**2 + w**2]
-    subs = _generate_substitutions(L, N,constraint, coords, size_tuple)
-    target_subs = [(y,-1), (x, -1-y**2)]
+    subs = _generate_substitutions(L, N, constraint, coords, size_tuple)
+    target_subs = [(y, -1), (x, -1 - y**2)]
 
     assert subs == target_subs
 
@@ -316,7 +315,7 @@ def test_cv_subs_func():
     bg = bgt.new()
     bg.add([c, se, kcl, r])
 
-    connect(c,(kcl,kcl.non_inverting))
+    connect(c, (kcl, kcl.non_inverting))
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
@@ -325,7 +324,7 @@ def test_cv_subs_func():
     subs = [(sympy.Symbol('u_0'), sympy.sympify('-exp(-t)'))]
 
     mappings, coords = inverse_coord_maps(*bg.basis_vectors)
-    assert _generate_cv_substitutions(cv_s, mappings,coords) == subs
+    assert _generate_cv_substitutions(cv_s, mappings, coords) == subs
 
 
 def test_cv_subs_const():
@@ -336,7 +335,7 @@ def test_cv_subs_const():
     bg = bgt.new()
     bg.add([c, se, kcl, r])
 
-    connect(c,(kcl,kcl.non_inverting))
+    connect(c, (kcl, kcl.non_inverting))
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
@@ -345,7 +344,7 @@ def test_cv_subs_const():
     subs = [(sympy.Symbol('u_0'), sympy.S(2))]
 
     mappings, coords = inverse_coord_maps(*bg.basis_vectors)
-    assert _generate_cv_substitutions(cv_s, mappings,coords) == subs
+    assert _generate_cv_substitutions(cv_s, mappings, coords) == subs
 
 
 def test_cv_subs_state_func():
@@ -356,7 +355,7 @@ def test_cv_subs_state_func():
     bg = bgt.new()
     bg.add([c, se, kcl, r])
 
-    connect(c,(kcl,kcl.non_inverting))
+    connect(c, (kcl, kcl.non_inverting))
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
@@ -365,8 +364,7 @@ def test_cv_subs_state_func():
     subs = [(sympy.Symbol('u_0'), sympy.sympify('-exp(-x_0)'))]
 
     mappings, coords = inverse_coord_maps(*bg.basis_vectors)
-    assert _generate_cv_substitutions(cv_s, mappings,coords) == subs
-
+    assert _generate_cv_substitutions(cv_s, mappings, coords) == subs
 
 
 def test_ported_cr():
@@ -393,13 +391,14 @@ def test_ported_cr():
     assert len(coords) == 15
 
     coords, mappings, lin_op, nl_op, conttr = model.system_model()
-    assert nl_op.is_zero
+    assert nl_op.is_zero_matrix
     assert not conttr
 
     assert model.constitutive_relations == [
         sympy.sympify('e_0 - 2*f_0 - 2*u_0')
     ]
 
+
 def test_ported_series_resistor():
 
     Se = new("Se")
@@ -409,7 +408,7 @@ def test_ported_series_resistor():
     ss = new("SS")
     model = new()
     model.add(
-        Se,r1,r2,kvl, ss
+        Se, r1, r2, kvl, ss
     )
     expose(ss)
     connect(Se, kvl.non_inverting)
@@ -423,6 +422,7 @@ def test_ported_series_resistor():
         sympy.sympify("e_0 - 3*f_0 - u_0")
     ]
 
+
 def test_ported_cap():
     model = new()
     c = new("C", value=3)
@@ -438,10 +438,11 @@ def test_ported_cap():
     expose(ss)
     assert len(model.ports) == 1
 
-
     assert model.constitutive_relations == [
-        sympy.sympify("dx_0 - f_0"),sympy.sympify("e_0 - x_0/3")
+        sympy.sympify("dx_0 - f_0"), sympy.sympify("e_0 - x_0/3")
     ]
+
+
 def test_ported_parallel_rc():
 
     model = new()
@@ -450,11 +451,11 @@ def test_ported_parallel_rc():
     zero = new("0")
     ss = new("SS")
     model.add(
-        r,c,zero, ss
+        r, c, zero, ss
     )
 
-    connect(r,zero)
-    connect(c,zero)
+    connect(r, zero)
+    connect(c, zero)
     connect(ss, zero)
 
     expose(ss)
@@ -465,6 +466,7 @@ def test_ported_parallel_rc():
         sympy.sympify("e_0 - x_0/3")
     ]
 
+
 class TestConstitutiveRelations:
 
     def test_empty_failstate(self):
diff --git a/test/test_biochem.py b/test/test_biochem.py
index d73df50..47f9c07 100644
--- a/test/test_biochem.py
+++ b/test/test_biochem.py
@@ -5,7 +5,7 @@
 import BondGraphTools as bgt
 from BondGraphTools import connect
 from BondGraphTools.exceptions import InvalidPortException
-from BondGraphTools.algebra import extract_coefficients, inverse_coord_maps,get_relations_iterator
+from BondGraphTools.algebra import extract_coefficients, inverse_coord_maps, get_relations_iterator
 from BondGraphTools.reaction_builder import Reaction_Network
 
 import logging
@@ -61,9 +61,10 @@ def test_re_con_rel():
 
     for r in get_relations_iterator(Re, mappings, coords):
         assert r in [
-            ({1:1, 3:1}, 0), ({1:1}, sympy.sympify("-r*exp(e_0) + r*exp(e_1)"))
+            ({1: 1, 3: 1}, 0), ({1: 1}, sympy.sympify("-r*exp(e_0) + r*exp(e_1)"))
         ]
 
+
 def test_a_to_b_model():
     A = bgt.new("Ce", library="BioChem", value=[1, 1, 1])
     B = bgt.new("Ce", library="BioChem", value=[1, 1, 1])
@@ -74,19 +75,20 @@ def test_a_to_b_model():
     Y_B = bgt.new('1')
 
     a_to_b = bgt.new()
-    a_to_b.add(A , Re, B,Y_A, Y_B)
+    a_to_b.add(A, Re, B, Y_A, Y_B)
 
     connect(A, Y_A.non_inverting)
     connect(B, Y_B.non_inverting)
-    connect((Re,0), Y_A.inverting)
-    connect((Re,1), Y_B.inverting)
+    connect((Re, 0), Y_A.inverting)
+    connect((Re, 1), Y_B.inverting)
 
-    eqns ={
+    eqns = {
         sympy.sympify("dx_0 + x_0 -x_1"), sympy.sympify("dx_1 + x_1 -x_0")
     }
     for relation in a_to_b.constitutive_relations:
         assert relation in eqns
 
+
 def test_ab_to_c_model():
 
     A = bgt.new("Ce", library="BioChem", value=[1, 1, 1])
@@ -104,14 +106,14 @@ def test_ab_to_c_model():
     connect(Re, Y_C)
     connect(Y_C, C)
 
-    state_basis ,_,_ = bg.basis_vectors
+    state_basis, _, _ = bg.basis_vectors
     (x_0, dx_0), = (k for k, (v, _) in state_basis.items() if v is A)
     (x_1, dx_1), = (k for k, (v, _) in state_basis.items() if v is B)
     (x_2, dx_2), = (k for k, (v, _) in state_basis.items() if v is C)
 
-    eqns = {dx_0 + x_0*x_1 - x_2,
-            dx_1 + x_0*x_1 - x_2,
-            dx_2 - x_0*x_1 + x_2}
+    eqns = {dx_0 + x_0 * x_1 - x_2,
+            dx_1 + x_0 * x_1 - x_2,
+            dx_2 - x_0 * x_1 + x_2}
 
     relations = set(bg.constitutive_relations)
     assert not relations ^ eqns
@@ -123,10 +125,10 @@ def test_new_reaction_network():
     rn.add_reaction(
         "A+B=C"
     )
-    assert rn.species ==["A","B","C"]
-    assert rn.forward_stoichiometry == sympy.Matrix([[1],[1],[0]])
+    assert rn.species == ["A", "B", "C"]
+    assert rn.forward_stoichiometry == sympy.Matrix([[1], [1], [0]])
     assert rn.reverse_stoichiometry == sympy.Matrix([[0], [0], [1]])
-    assert rn.stoichiometry == sympy.Matrix([[-1],[-1],[1]])
+    assert rn.stoichiometry == sympy.Matrix([[-1], [-1], [1]])
 
 
 def test_rn_to_bond_graph():
@@ -146,7 +148,7 @@ def test_cat_rn():
     assert len(rn._reactions) == 2
 
 
-#TODO: fix when we rework parameters
+# TODO: fix when we rework parameters
 def test_nlin_se():
     rn = Reaction_Network(name="A+B to C", reactions="A+B=C")
     system = rn.as_network_model(normalised=True)
@@ -154,7 +156,7 @@ def test_nlin_se():
     for param in system.params:
         system.set_param(param, 1)
 
-    Ce_A = system/ "A"
+    Ce_A = system / "A"
     Ce_B = system / "B"
     Ce_C = system / "C"
 
@@ -174,11 +176,11 @@ def test_nlin_se():
     connect(Se_A, J_A)
 
     assert len(system.state_vars) == 3
-    state_basis , _, _ = system.basis_vectors
+    state_basis, _, _ = system.basis_vectors
 
-    (x0, dx0), = (k for k, (v,_) in  state_basis.items() if v is Ce_A)
-    (x1, dx1), = (k for k, (v,_) in  state_basis.items() if v is Ce_B)
-    (x2, dx2), = (k for k, (v,_) in  state_basis.items() if v is Ce_C)
+    (x0, dx0), = (k for k, (v, _) in state_basis.items() if v is Ce_A)
+    (x1, dx1), = (k for k, (v, _) in state_basis.items() if v is Ce_B)
+    (x2, dx2), = (k for k, (v, _) in state_basis.items() if v is Ce_C)
     E = sympy.S("E")
     solutions = {dx1 + E * x1 - x2, dx2 - E * x1 + x2, x0 - E, dx0}
 
@@ -202,5 +204,3 @@ def test_naming(self):
 
         assert isinstance(re, SymmetricComponent)
         assert re.metamodel == 'R'
-
-
diff --git a/test/test_component_manager.py b/test/test_component_manager.py
index aadd29e..06241eb 100644
--- a/test/test_component_manager.py
+++ b/test/test_component_manager.py
@@ -6,7 +6,7 @@
 
 def test_load():
 
-    dir, _ = os.path.split(__file__ )
+    dir, _ = os.path.split(__file__)
     file = os.path.join(dir, "lib/test_lib.json")
 
     assert cm.load_library(file)
@@ -16,7 +16,7 @@ def test_load():
 def test_duplicate_load():
 
     del cm.__libraries["test"]
-    dir, _ = os.path.split(__file__ )
+    dir, _ = os.path.split(__file__)
     file = os.path.join(dir, "lib/test_lib.json")
 
     assert cm.load_library(file)
@@ -30,5 +30,3 @@ def test_bad_libary():
     file = os.path.join(dir, "lib/test_lib2.json")
 
     assert not cm.load_library(file)
-
-
diff --git a/test/test_core_functionality.py b/test/test_core_functionality.py
index 26d4c28..9a39795 100644
--- a/test/test_core_functionality.py
+++ b/test/test_core_functionality.py
@@ -37,7 +37,6 @@ def test_add():
     bg = bgt.new()
     bg.add([c, se])
 
-
     assert len(bg.internal_ports) == 2
     assert len(bg.state_vars) == 1
 
@@ -58,7 +57,6 @@ def test_equal():
     assert c_1 is not c_2
 
 
-
 class TestBond:
     def test_create(self):
         from BondGraphTools.base import Bond
@@ -67,7 +65,7 @@ def test_create(self):
         c = new('C')
         one = new('1')
 
-        b_1 = Bond(head=(c,0), tail=(one, 0))
+        b_1 = Bond(head=(c, 0), tail=(one, 0))
 
         assert c in b_1
         assert one in b_1
@@ -124,13 +122,13 @@ def test_one_port(self):
         r = new('R')
         one = new('1')
         bg = new()
-        bg.add(c,se,r,one)
+        bg.add(c, se, r, one)
 
         connect(c, one.non_inverting)
         connect(se, one)
         connect(one, r)
 
-        p0,p1,p2 = tuple(one.ports)
+        p0, p1, p2 = tuple(one.ports)
 
         assert p0.weight == 1
         assert p1.weight == 1
@@ -167,10 +165,11 @@ def test_component_in_bond(self):
                  (se, one),
                  (one, r)]
 
-        for bond in bonds: connect(*bond)
+        for bond in bonds:
+            connect(*bond)
 
-        comps = [{c, one}, {se, one}, {one ,r}]
-        all_comps = {c,se,r,one}
+        comps = [{c, one}, {se, one}, {one, r}]
+        all_comps = {c, se, r, one}
 
         for i, bond in enumerate(bg.bonds):
             for component in all_comps:
@@ -179,8 +178,6 @@ def test_component_in_bond(self):
                 else:
                     assert component not in bond
 
-
-
     def test_component_in_bond(self):
         # see issues 85
         c = new('C')
@@ -194,10 +191,11 @@ def test_component_in_bond(self):
                  (se, one),
                  (one, r)]
 
-        for bond in bonds: connect(*bond)
+        for bond in bonds:
+            connect(*bond)
 
-        comps = [{c, one}, {se, one}, {one ,r}]
-        all_comps = {c,se,r,one}
+        comps = [{c, one}, {se, one}, {one, r}]
+        all_comps = {c, se, r, one}
 
         for i, bond in enumerate(bg.bonds):
             for component in all_comps:
@@ -206,6 +204,7 @@ def test_component_in_bond(self):
                 else:
                     assert component not in bond
 
+
 def test_disconnect_ports():
 
     c = bgt.new("C")
@@ -226,6 +225,7 @@ def test_disconnect_ports():
     disconnect(c, se)
     assert not bg.bonds
 
+
 def test_many_port():
 
     j = bgt.new("0")
@@ -257,7 +257,7 @@ def test_delete_one_from_many_port():
     connect(c, j)
     connect(se, j)
 
-    (p,q), (r,s) = bg.bonds[1]
+    (p, q), (r, s) = bg.bonds[1]
 
     assert p in (se, j)
     assert r in (se, j)
@@ -319,13 +319,13 @@ def test_disconnect_multiport():
     bg.add([zero, r, c, one])
 
     connect(zero, one.non_inverting)
-    connect(r,zero)
+    connect(r, zero)
     connect(c, one.inverting)
 
     assert len(bg.bonds) == 3
     disconnect(zero, one)
     assert len(bg.bonds) == 2
-    assert ((zero, 0), (one,0)) not in bg.bonds
+    assert ((zero, 0), (one, 0)) not in bg.bonds
 
 
 def test_disconnect_component():
@@ -360,7 +360,7 @@ def test_remove_component(self):
         connect(c, zero)
         r_p, = r.ports
         c_p, = c.ports
-        z0,z1, = zero.ports
+        z0, z1, = zero.ports
         assert c in bg.components
         assert set(bg.bonds) == {
             (r_p, z0),
@@ -391,9 +391,9 @@ def test_swap_component_1(self):
         c = bgt.new("C", value=1)
 
         bg = bgt.new()
-        bg.add([zero, r, c ])
+        bg.add([zero, r, c])
 
-        connect(r,zero)
+        connect(r, zero)
         connect(c, zero)
         assert c in bg.components
         r_p, = r.ports
@@ -413,7 +413,6 @@ def test_swap_component_1(self):
         assert len(bg.state_vars) == 0
         assert len(bg.control_vars) == 1
 
-
         assert bg.bonds == [
             (r_p, z0),
             (sf_port, z1)
@@ -424,7 +423,6 @@ def test_swap_component_1(self):
 
         swap(Sf, c)
 
-
         assert bg.bonds == [
             (r_p, z0),
             (c_p, z1)
@@ -440,9 +438,9 @@ def test_swap_components_2(self):
         c = bgt.new("C")
 
         bg = bgt.new()
-        bg.add([zero, r, c ])
+        bg.add([zero, r, c])
 
-        connect(r,zero)
+        connect(r, zero)
         connect(c, zero)
 
         p1, = r.ports
@@ -463,7 +461,7 @@ def test_swap_failstate(self):
         l = bgt.new("I")
 
         bg = bgt.new()
-        bg.add([zero, r, c ])
+        bg.add([zero, r, c])
         connect(r, zero)
         connect(c, zero)
 
diff --git a/test/test_diode.py b/test/test_diode.py
index 239582c..3e01372 100644
--- a/test/test_diode.py
+++ b/test/test_diode.py
@@ -6,11 +6,11 @@
 
 def test_make_diode():
 
-    d = bgt.new("Di", library="elec", value={"Is": 5*10**-6})
+    d = bgt.new("Di", library="elec", value={"Is": 5 * 10**-6})
     # todo: fix this
     # assert d.metamodel == "Di"
     assert len(d.ports) == 1
-    assert d.params["Is"]["value"] == 5*10**-6
+    assert d.params["Is"]["value"] == 5 * 10**-6
 
 
 def test_diode_model():
@@ -21,4 +21,3 @@ def test_diode_model():
     assert not d.control_vars
 
     assert d.constitutive_relations[0] == sp.sympify('f_0 - exp(e_0) +1')
-
diff --git a/test/test_io.py b/test/test_io.py
index 809640f..fb0804c 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -10,6 +10,7 @@
 file_path = pathlib.Path(__file__).parent / 'files'
 logging.basicConfig(level=logging.DEBUG)
 
+
 def test_load_rlc():
 
     path = str(file_path / 'rlc.bg')
@@ -21,21 +22,21 @@ def test_load_rlc():
             "RLC:/kcl",
             "RLC:/Sf"]
 
-    c, r, l, kcl, sf =  (comp for uri in uris for comp in model.components if
-                      comp.uri == uri)
+    c, r, l, kcl, sf = (comp for uri in uris for comp in model.components if
+                        comp.uri == uri)
 
     assert len(model.control_vars) == 1
     assert 'u_0' in model.control_vars
 
     assert c.params['C']['value'] == 10
-    assert r.params['r']['value']  == 100
-    assert l.params['L']['value']  == 10
+    assert r.params['r']['value'] == 100
+    assert l.params['L']['value'] == 10
 
     r_0, = r.ports
     c_0, = c.ports
     l_0, = l.ports
-    sf_0,  = sf.ports
-    kcl_0, kcl_1,kcl_2,kcl_3 = kcl.ports
+    sf_0, = sf.ports
+    kcl_0, kcl_1, kcl_2, kcl_3 = kcl.ports
 
     assert set(model.bonds) == {
         (r_0, kcl_0),
@@ -55,7 +56,7 @@ def test_load_rlc():
     }
 
     assert (set(model.constitutive_relations) == eqns) or \
-                (set(model.constitutive_relations) == eqns_2)
+        (set(model.constitutive_relations) == eqns_2)
 
 
 def test_load_rlc_parallel():
@@ -81,6 +82,7 @@ def test_load_rlc_parallel():
     for r in rel:
         assert r in (eq1, eq2)
 
+
 def test_load_modular():
 
     model_1 = load(file_path / "modular.bg")
@@ -89,6 +91,7 @@ def test_load_modular():
     assert model_1.name == "system"
 
     tree = set()
+
     def uri_tree(bg):
         tree.add(bg.uri)
         try:
@@ -126,7 +129,7 @@ def test_load_model_from_modular():
     assert len(Vs.components) == 3
     assert len(Vs2.components) == 3
 
-@pytest.mark.usefixture("rlc")
+
 def test_save_build_component(rlc):
 
     r, = (c for c in rlc.components if c.metamodel == "R")
@@ -139,7 +142,7 @@ def test_save_build_component(rlc):
     assert r_str == f"{r.name} base/R r=10"
     assert c_str == f"{c.name} base/C"
 
-@pytest.mark.usefixture("rlc")
+
 def test_save_build_model(rlc):
     r, = (c for c in rlc.components if c.metamodel == "R")
     c, = (c for c in rlc.components if c.metamodel == "C")
@@ -158,6 +161,7 @@ def test_save_build_model(rlc):
         f"{r.name} {kvl.name}", f"{c.name} {kvl.name}", f"{l.name} {kvl.name}"
     }
 
+
 def test_build_templated_model():
     root = new(name="System")
     model = new(name="Vs")
@@ -165,8 +169,8 @@ def test_build_templated_model():
     ss = new("SS", name="pout")
     Se = new("Se", name="Vs")
     zero = new("0", name="0_0")
-    model.add(ss,Se,zero)
-    connect(ss,zero)
+    model.add(ss, Se, zero)
+    connect(ss, zero)
     connect(Se, zero)
 
     file = file_path / "temp.bg"
@@ -177,7 +181,7 @@ def test_build_templated_model():
         save(root, file)
 
         with open(file, 'r') as fs:
-            temp_data = yaml.load(fs)
+            temp_data = yaml.load(fs, Loader=yaml.SafeLoader)
 
     assert temp_data["root"] == root.name
     assert temp_data["models"].keys() == {
@@ -193,7 +197,6 @@ def test_build_templated_model():
     }
 
 
-@pytest.mark.usefixture("rlc")
 def test_rlc_save(rlc):
     filename = str(file_path / "test_rlc.bg")
     rlc.name = "RLC"
@@ -202,15 +205,15 @@ def test_rlc_save(rlc):
     l, = (c for c in rlc.components if c.metamodel == "I")
     kvl, =  (c for c in rlc.components if c.metamodel == "0")
 
-    r.name= "R1"
-    c.name= "C1"
-    l.name= "L1"
+    r.name = "R1"
+    c.name = "C1"
+    l.name = "L1"
 
     with TempFile(filename):
         save(rlc, filename)
 
         with open(filename, 'r') as fs:
-            test_data = yaml.load(fs)
+            test_data = yaml.load(fs, Loader=yaml.SafeLoader)
 
     assert test_data["version"] == dm.FILE_VERSION
     assert test_data["root"] == "RLC"
@@ -220,7 +223,7 @@ def test_rlc_save(rlc):
     assert set(model_data.keys()) == {"components", "netlist"}
 
     assert set(model_data["components"]) == {
-        "R1 base/R r=1", "C1 base/C C=1", "L1 base/I L=1","kvl base/0"}
+        "R1 base/R r=1", "C1 base/C C=1", "L1 base/I L=1", "kvl base/0"}
 
     assert set(model_data["netlist"]) == {
         "R1 kvl", "C1 kvl", "L1 kvl"
@@ -239,8 +242,11 @@ def test_save_idempotent():
     assert model_1.name == model_2.name
     assert_components_are_equal(model_1, model_2)
 
+
 def assert_components_are_equal(bg_1, bg_2):
-    names = {c.name for c in bg_1.components} & {c.name for c in bg_2.components}
+    names = {
+        c.name for c in bg_1.components} & {
+        c.name for c in bg_2.components}
     if len(names) != len(bg_1.components):
         assert False, names
 
@@ -250,9 +256,10 @@ def assert_components_are_equal(bg_1, bg_2):
         if hasattr(c1, 'components'):
             assert_components_are_equal(c1, c2)
         else:
-            assert_atomics_are_equal(c1,c2)
+            assert_atomics_are_equal(c1, c2)
 
-def assert_atomics_are_equal(c1,c2):
+
+def assert_atomics_are_equal(c1, c2):
     ports_1 = {repr(p) for p in c1.ports}
     ports_2 = {repr(p) for p in c2.ports}
     assert ports_1 == ports_2
@@ -276,9 +283,9 @@ def __enter__(self):
             pass
 
     def __exit__(self, exc_type, exc_val, exc_tb):
-        with open(self.file, 'r' ) as file:
+        with open(self.file, 'r') as file:
             print("".join(line for line in file.readlines()))
         try:
             os.remove(self.file)
         except FileNotFoundError:
-            pass
\ No newline at end of file
+            pass
diff --git a/test/test_layout.py b/test/test_layout.py
index 8648636..0c95a1c 100644
--- a/test/test_layout.py
+++ b/test/test_layout.py
@@ -2,11 +2,8 @@
 from BondGraphTools.view import _build_graph
 
 
-@pytest.mark.usefixture("rlc")
 def test_edge_list(rlc):
 
     graph = _build_graph(rlc)
 
     assert graph.nnz == 6
-
-
diff --git a/test/test_modular.py b/test/test_modular.py
index b71424a..bdbf0f1 100644
--- a/test/test_modular.py
+++ b/test/test_modular.py
@@ -24,13 +24,14 @@ def test_source_sensor():
     expose(ss, "A")
     assert len(model.ports) == 1
 
+
 def test_ss_exposure():
     model = new()
     ss = new("SS")
     sf = new("Sf", name="Sf")
     zero = new("0")
 
-    model.add(ss,sf,zero)
+    model.add(ss, sf, zero)
 
     connect(sf, zero)
     connect(ss, zero)
@@ -44,12 +45,13 @@ def test_ss_exposure():
     assert p.name is "pin"
     assert set(model.control_vars.values()) == {(sf, 'f')}
 
+
 def test_load_modular():
     model_1 = load(file_path / "modular.bg")
 
     assert model_1.name == "system"
 
-    assert {c.name for c in model_1.components}  == {"Vs", "Z", "kvl"}
+    assert {c.name for c in model_1.components} == {"Vs", "Z", "kvl"}
     for c in model_1.components:
         assert c
         assert c.parent is model_1
@@ -67,7 +69,6 @@ def test_load_modular():
     assert len(Vs.control_vars) == 1
 
 
-
 def test_modularity():
     model_1 = load(file_path / "modular.bg")
 
@@ -75,14 +76,10 @@ def test_modularity():
     Z = model_1 / "Z"
 
     Vs_cr = Vs.constitutive_relations
-    assert Vs_cr ==[sp.sympify('f_0 + u_0')]
-    
+    assert Vs_cr == [sp.sympify('f_0 + u_0')]
+
     assert (set(model_1.constitutive_relations) == {
         sp.sympify("dx_0 - x_1"), sp.sympify("dx_1 - u_0 +x_0 + x_1")
-    }) or ( set(model_1.constitutive_relations) == {
+    }) or (set(model_1.constitutive_relations) == {
         sp.sympify("dx_1 - x_0"), sp.sympify("dx_0 - u_0 +x_0 + x_1")}
-           )
-
-
-
-
+    )
diff --git a/test/test_porthamiltonian.py b/test/test_porthamiltonian.py
index 50766c6..4ddcccc 100644
--- a/test/test_porthamiltonian.py
+++ b/test/test_porthamiltonian.py
@@ -11,15 +11,15 @@ def test_hamiltonian():
     )
 
     assert state_vars == {"q_0": "x"}
-    assert params == {"C": None, "D":None}
-    assert ports == {0:None}
+    assert params == {"C": None, "D": None}
+    assert ports == {0: None}
     assert relations == ["-e_0 + q_0**3/D + q_0/C", "dq_0 - f_0"]
 
 
 def test_create_PH():
     hamiltonian = "x^2/(2*C) + x^4/(4*D)"
     ph = bgt.new(component="PH", value=hamiltonian)
-    assert ph.params == {"C": None, "D":None}
+    assert ph.params == {"C": None, "D": None}
 
     assert ph.state_vars == {"q_0": "x"}
 
@@ -49,7 +49,6 @@ def test_create_PH_2():
                 "dq_1 - f_1"):
         assert sp.sympify(rel) in ph.constitutive_relations
 
-
     assert ph.hamiltonian == hamiltonian
 
     port_list = list(ph.ports)
@@ -64,7 +63,7 @@ def test_create_PH_parameters():
     p_2 = sp.S("k")
     build_args = {
         "hamiltonian": hamiltonian,
-        "params":{"C":p_1, "D":p_2}
+        "params": {"C": p_1, "D": p_2}
     }
 
     ph = bgt.new(component="PH", value=build_args)
diff --git a/test/test_ports.py b/test/test_ports.py
index 6c67e47..980a84b 100644
--- a/test/test_ports.py
+++ b/test/test_ports.py
@@ -1,5 +1,6 @@
 import pytest
 
+
 class TestPort:
     def test_in(self):
         from BondGraphTools.base import Port
@@ -28,6 +29,7 @@ def test_cmp(self):
         assert port_c is not port_one
         assert port_c != port_one
 
+
 class TestIntegrations:
 
     def test_get_exposed_port(self):
@@ -53,10 +55,10 @@ def test_failstate(self):
         expose(port, label='port_1')
 
         test_model = new()
-        component = new("0") ## could be anything
+        component = new("0")  # could be anything
         test_model.add(component, model)
 
-        target_port =(model, 'port_')     # typo in port name
+        target_port = (model, 'port_')     # typo in port name
         with pytest.raises(InvalidPortException):
             connect(component, target_port)
 
@@ -76,4 +78,3 @@ def test_connection_feedback(self):
             connect((c2, 0), c3)
 
         assert "Port is already connected: " in str(ex)
-
diff --git a/test/test_rtd_examples.py b/test/test_rtd_examples.py
index e4571cc..ea9da1d 100644
--- a/test/test_rtd_examples.py
+++ b/test/test_rtd_examples.py
@@ -23,10 +23,19 @@ def enzyme_catalysed_reaction(name):
         # Here we build the reaction, again with the rate as a control variable.
         # Again, we assume parameterised have be normalised with respect to
         # pressure and temperature.
-        cat_reaction = bgt.new("Re", name="Re", library="BioChem", value={'r': None, 'R': 1, 'T': 1})
+        cat_reaction = bgt.new(
+            "Re", name="Re", library="BioChem", value={
+                'r': None, 'R': 1, 'T': 1})
 
         # We choose 'k' to be 1 for demonstration.
-        enzyme = bgt.new("Ce", name="E", library="BioChem", value={'k': 1, 'R': 1, 'T': 1})
+        enzyme = bgt.new(
+            "Ce",
+            name="E",
+            library="BioChem",
+            value={
+                'k': 1,
+                'R': 1,
+                'T': 1})
 
         # Substrate + Enzyme flux conservation law
         SE = bgt.new('1')
@@ -36,7 +45,15 @@ def enzyme_catalysed_reaction(name):
         # Conservation of enzyme law.
         law_E = bgt.new("0")
 
-        bgt.add(cat_model, substrate, product, enzyme, SE, PE, law_E, cat_reaction)
+        bgt.add(
+            cat_model,
+            substrate,
+            product,
+            enzyme,
+            SE,
+            PE,
+            law_E,
+            cat_reaction)
 
         connections = [
             (substrate, SE),
@@ -55,4 +72,4 @@ def enzyme_catalysed_reaction(name):
 
         return cat_model
 
-    E1 = enzyme_catalysed_reaction("E1")
\ No newline at end of file
+    E1 = enzyme_catalysed_reaction("E1")
diff --git a/test/test_sims.py b/test/test_sims.py
index 06f7dc7..608a8eb 100644
--- a/test/test_sims.py
+++ b/test/test_sims.py
@@ -2,19 +2,18 @@
 import numpy as np
 import sympy as sp
 from BondGraphTools import *
-from BondGraphTools.config import config
 from BondGraphTools.exceptions import ModelException
-from BondGraphTools.sim_tools import simulate, to_julia_function_string
+from BondGraphTools.sim_tools import simulate, _bondgraph_to_residuals
 from BondGraphTools.algebra import inverse_coord_maps
 
-@pytest.mark.slow
+
 def test_c_sim_fail():
 
     c = new("C")
     with pytest.raises(ModelException):
-        t, x = simulate(c, timespan=[0, 1], x0=[1],dx0=[1])
+        t, x = simulate(c, timespan=[0, 1], x0=[1], dx0=[1])
+
 
-@pytest.mark.slow
 def test_c_se_build_ode():
 
     c = new("C", value=1)
@@ -24,21 +23,33 @@ def test_c_se_build_ode():
     bg = new()
     bg.add([c, se, kcl, r])
 
-    connect(c,(kcl, kcl.non_inverting))
+    connect(c, (kcl, kcl.non_inverting))
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
     # "dx_0 - u_0 + x_0"
     # so f(x,t) = exp(-t) - x
 
-    j = config.julia
-    func_str, diff_vars = to_julia_function_string(bg,
-                                                   control_vars=['-exp(-t)'],
-                                                   in_place=False)
-    func = j.eval(func_str)
-    assert func(0, 0, 0, 0) == -1
-    assert func(2, 0, 0, 0) == 1
-    assert func(0, 2, 0, 0) == 1
+    def u(t, x, dx):
+        return -np.exp(-t)
+
+    residual_func, diff_vars = _bondgraph_to_residuals(bg, control_vars=[u])
+
+    r = [0]
+    residual_func(0, [0], [0], r)
+    assert r == [-1]
+    residual_func(0, [0], [2], r)
+    assert r == [1]
+    residual_func(0, [2], [0], r)
+    assert r == [1]
+
+    residual_func(0, [0], [1], r)
+    assert r == [0]
+
+    t_test = np.log(4)
+    residual_func(t_test, [1 / 8], [1 / 8], r)
+    assert r == [0]
+
 
 @pytest.mark.slow
 def test_c_se_sim():
@@ -54,26 +65,28 @@ def test_c_se_sim():
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
+    def u(t, x, dx):
+        return -np.exp(-t)
+    assert str(bg.constitutive_relations) == '[dx_0 + u_0 + x_0]'
     with pytest.raises(ModelException) as ex:
         t, x = simulate(
-            c, timespan=[0, 10], x0=[0], dx0=[1]
+            c, timespan=[0, 10], x0=[0]
         )
-        assert "Control variable not specified" in ex.args
+        assert "Control variable u_0 must be specified" in ex.args
     t, x = simulate(
-        bg, timespan=[0, 10], x0=[0], dx0=[1], control_vars=['-exp(-t)']
+        bg, timespan=[0, 10], x0=[0], control_vars=[u]
     )
 
     assert t[0] == 0
     assert t[-1] == 10
-    t_cols, _ = t.shape
-    assert t_cols > 1
-    assert (t_cols, 1) == x.shape
+
+    assert (len(t), 1) == x.shape
     assert x[0, 0] == 0
+    solution = t * np.exp(-t)
 
-    solution = (1+t)*np.exp(-t)
-    solution = solution.reshape(t_cols, 1)
+    res = abs(x - solution)
+    assert (res < 0.001).all()
 
-    assert ((x - solution) < 0.001).all()
 
 @pytest.mark.slow
 def test_c_se_sum_switch():
@@ -88,14 +101,16 @@ def test_c_se_sum_switch():
     connect(r, (kcl, kcl.non_inverting))
     connect(se, (kcl, kcl.non_inverting))
 
-    bang_bang = ["x_0 >= 1 ?  1.5: -2 "]
+    def bang_bang(t, x, dx):
+        return 1.5 if x >= 1 else -2.0
 
     t, x = simulate(
-        bg, timespan=[0, 10], x0=[0],dx0=[1], control_vars=bang_bang
+        bg, timespan=[0, 10], x0=[0], dx0=[1], control_vars=[bang_bang]
     )
 
     assert (x[0, -1] - 1) < 0.001
 
+
 @pytest.mark.skip
 @pytest.mark.slow
 def test_rlc():
@@ -108,13 +123,11 @@ def test_rlc():
     bg = new()
     bg.add([c, se, kvl, r, l])
 
-
     connect(c, kvl)
     connect(r, kvl)
     connect(se, kvl)
     connect(l, kvl)
 
     t, x = simulate(
-        bg, timespan=[0, 10], x0=[1,0], control_vars=[1]
+        bg, timespan=[0, 10], x0=[1, 0], control_vars=[1]
     )
-
diff --git a/test/test_uri.py b/test/test_uri.py
index 0ba6b26..4e596ee 100644
--- a/test/test_uri.py
+++ b/test/test_uri.py
@@ -24,15 +24,12 @@ def test_uri_tree():
 
     assert c_o.uri == "Model:/C_1"
 
-
     assert internal_model in model.components
     assert r_i in internal_model.components
     assert r_i not in model.components
 
-
     assert r_i.uri == "Model:/Model/" + r_i.name
 
-
     c_i = new("C", name="C_1")
     internal_model.add(
         c_i
@@ -47,4 +44,3 @@ def test_uri_tree():
 
     assert c_i.uri == "Model:/C_1"
     assert c_o.uri == "Model:/C_1"
-