diff --git a/README.md b/README.md
index 86a4bc7..e97867b 100644
--- a/README.md
+++ b/README.md
@@ -43,8 +43,9 @@ Systems I have tested with this code:
 - Lorenz Attractor dynamical system
 
 To run this module, you need a system of 3 dynamical ODE's.  This module does *not* use SymPy, so make sure to
-verify that you've typed in the equations correctly.  Input the 3 functions as f1, f2, and f3, and then call 
-`abm_comp` for a 3D plot of the solution.
+verify that you've typed in the equations correctly.  Input the 3 functions as f1, f2, and f3, and then
+instantiate the class `ABM` as shown below.  To 3D plot. use ABM.plot3d().  To get numerical answers,
+use ABM.adams().
 
 This is the Rossler Attractor:
 
@@ -58,7 +59,7 @@ def f2(x, y, z):
 def f3(x, y, z):
     return 0.1 + z*(x - 14)
 
-abm_comp(f1, f2, f3, 0, 15, 15, 36, 0, 100, 10000)
+res = ABM(f1, f2, f3, 0, 15, 15, 36, 0, 100, 10000)
 # func1, func2, func3, initial t0, x0, y0, z0, lower bound, upper, num iterations
 ```
 
diff --git a/computepac/__init__.py b/computepac/__init__.py
index b80e9a4..44b7b8b 100644
--- a/computepac/__init__.py
+++ b/computepac/__init__.py
@@ -6,5 +6,5 @@
 from .dynamic_triple import *
 from computepac import *
 
-__author__ = 'mathemacode'
-__copyright__ = 'License in LICENSE.txt'
+__author__ = "mathemacode"
+__copyright__ = "License in LICENSE.txt"
diff --git a/computepac/bisection.py b/computepac/bisection.py
index fa2be1c..3e529d9 100644
--- a/computepac/bisection.py
+++ b/computepac/bisection.py
@@ -18,14 +18,14 @@
 
 
 def bisection(f, a, b, acc) -> list:
-    x = sym.Symbol('x')
-    start_time = time.time()                    # start timer
-    error = 10                                  # before assignment in while loop
-    err = []                                    # to store errors
+    x = sym.Symbol("x")
+    start_time = time.time()  # start timer
+    error = 10  # before assignment in while loop
+    err = []  # to store errors
 
-    while error > acc:                          # if error still larger than desired accuracy
+    while error > acc:  # if error still larger than desired accuracy
 
-        m = (a + b) / 2                         # middle between bounds
+        m = (a + b) / 2  # middle between bounds
         ans = f.evalf(subs={x: m})
 
         if ans > 0:
@@ -39,8 +39,8 @@ def bisection(f, a, b, acc) -> list:
         err.append(error)
 
     end_time = time.time()
-    total_time = end_time-start_time
+    total_time = end_time - start_time
     num_itr = len(err)
     root = m
-    
+
     return [num_itr, root, total_time]
diff --git a/computepac/dynamic_triple.py b/computepac/dynamic_triple.py
index bd23e11..cd30f59 100644
--- a/computepac/dynamic_triple.py
+++ b/computepac/dynamic_triple.py
@@ -7,8 +7,12 @@
 # Input t, x, y, z initial conditions
 # Input a, b, n as lower bound, upper bound, num iterations
 #
-# Can run plot(a, b, n) to see results.  Can also run
-# rk() or adams() by themselves, they output a list of lists
+# This file is a Class structure so you must instantiate it first,
+# ex: my_abm_calculation = ABM( enter all params here ), then call the
+# other methods ex.: my_abm_calculation.plot3d()
+#
+# Can run ABM.plot() to see results.  Can also run
+# ABM.rk() or ABM.adams() by themselves, they output a list of lists
 # result = [[xvals], [yvals], [zvals]] for plotting or analysis
 #
 # Note line in plot() that specifies RK vs. ABM results.  You can
@@ -29,175 +33,244 @@
 #        def f3(x, y, z):
 #            return 0.1 + z*(x - 14)
 #
-#        abm_comp(f1, f2, f3, 0, 15, 15, 36, 0, 100, 10000)
+#        res = ABM(f1, f2, f3, 0, 15, 15, 36, 0, 100, 10000)
+#        res.plot()
 #
-# Output will be 3D matplotlib plot if you do this
+# Output will be 3D matplotlib plot if you do this.  Other methods like res.adams() will
+# show numerical results.
 #
 
 import numpy as np
 import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-
-
-def abm_comp(f1, f2, f3, t0, x0, y0, z0, a, b, n):
-    # function 1, func2, func3, initial t0, x0, y0, z0, lower bound, upper, num iterations
-
-    global params
-    params = [t0, x0, y0, z0]
-
-    plot(a, b, n)
-
-
-# derivative functions
-def x1p(x1, x2, x3):
-    return f1(x1, x2, x3)
-
-
-def x2p(x1, x2, x3):
-    return f2(x1, x2, x3)
-
-
-def x3p(x1, x2, x3):
-    return f3(x1, x2, x3)
-
-
-def adams(a, b, n) -> list:  # A-B-M computation, using RK4 to start
-
-    h = (b - a) / n
-    y1 = np.zeros(n)  # initialize y1, y2, y3 matrices (and y_predictors)
-    y2 = np.zeros(n)
-    y3 = np.zeros(n)
-
-    y_pred1 = np.zeros(n)
-    y_pred2 = np.zeros(n)
-    y_pred3 = np.zeros(n)
-
-    for E in (0, 1, 2, 3):
-        y1[E] = rk(a, b, n, E, 1)  # initial estimations to start A-B-M
-        y2[E] = rk(a, b, n, E, 2)
-        y3[E] = rk(a, b, n, E, 3)
-
-    for i in range(4, n):  # Equations change slightly, hard to refactor to make clean
-
-        # x1p
-        y_pred1[i] = y1[i - 1] + (h / 24 * (55 *
-                                            x1p(y1[i - 1], y2[i - 1], y3[i - 1]) - 59 *
-                                            x1p(y1[i - 2], y2[i - 2], y3[i - 2]) + 37 *
-                                            x1p(y1[i - 3], y2[i - 3], y3[i - 3]) - 9 * 0))
-
-        # x1p
-        y1[i] = y1[i - 1] + (h / 24 * (9 *
-                                       x1p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1]) + 19 *
-                                       x1p(y1[i - 1], y2[i - 1], y3[i - 1]) - 5 *
-                                       x1p(y1[i - 2], y2[i - 2], y3[i - 2]) +
-                                       x1p(y1[i - 3], y2[i - 3], y3[i - 3])))
-
-        # x2p
-        y_pred2[i] = y2[i - 1] + (h / 24 * (55 *
-                                            x2p(y1[i - 1], y2[i - 1], y3[i - 1]) - 59 *
-                                            x2p(y1[i - 2], y2[i - 2], y3[i - 2]) + 37 *
-                                            x2p(y1[i - 3], y2[i - 3], y3[i - 3]) - 9 * (15 * (-8) - 15)))
-
-        # x2p
-        y2[i] = y2[i - 1] + (h / 24 * (9 *
-                                       x2p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1]) + 19 *
-                                       x2p(y1[i - 1], y2[i - 1], y3[i - 1]) - 5 *
-                                       x2p(y1[i - 2], y2[i - 2], y3[i - 2]) +
-                                       x2p(y1[i - 3], y2[i - 3], y3[i - 3])))
-
-        # x3p
-        y_pred3[i] = y3[i - 1] + (h / 24 * (55 *
-                                            x3p(y1[i - 1], y2[i - 1], y3[i - 1]) - 59 *
-                                            x3p(y1[i - 2], y2[i - 2], y3[i - 2]) + 37 *
-                                            x3p(y1[i - 3], y2[i - 3], y3[i - 3]) - 9 * ((15 * 15) - ((8 / 3) * 36))))
-
-        # x3p
-        y3[i] = y3[i - 1] + (h / 24 * (9 *
-                                       x3p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1]) + 19 *
-                                       x3p(y1[i - 1], y2[i - 1], y3[i - 1]) - 5 *
-                                       x3p(y1[i - 2], y2[i - 2], y3[i - 2]) +
-                                       x3p(y1[i - 3], y2[i - 3], y3[i - 3])))
-
-    return [y1, y2, y3]
-
-
-def rk(a, b, n, whichy=9, whichx=9) -> list:  # Runge-Kutta 4 computation by itself
-    # a lower bound
-    # b upper bound
-    # n number of steps
-    # whichy indicator for which Y1, Y2, Y3 val is needed to be return by function
-    # whichx means x1, x2, x3 when calling for values Y1, Y2, Y3, etc...
-
-    h = (b - a) / n  # step size
-    p = 0  # flag for while loop
 
-    # intialize lists for plot
-    t = np.zeros(int(n + 1))
-    x1 = np.zeros(int(n + 1))
-    x2 = np.zeros(int(n + 1))
-    x3 = np.zeros(int(n + 1))
 
-    # define initial values
-    t[0] = params[0]
-    x1[0] = params[1]
-    x2[0] = params[2]
-    x3[0] = params[3]
-
-    # initialize slots for values of RK4
-    i = np.zeros(4)  # for x1
-    j = np.zeros(4)  # for x2
-    k = np.zeros(4)  # for x3
-
-    # iterate n times
-    while p < n:
-        i[0] = h * x1p(x1[p], x2[p], x3[p])
-        j[0] = h * x2p(x1[p], x2[p], x3[p])
-        k[0] = h * x3p(x1[p], x2[p], x3[p])
-
-        i[1] = h * x1p(x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0])
-        j[1] = h * x2p(x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0])
-        k[1] = h * x3p(x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0])
-
-        i[2] = h * x1p(x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1])
-        j[2] = h * x2p(x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1])
-        k[2] = h * x3p(x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1])
-
-        i[3] = h * x1p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
-        j[3] = h * x2p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
-        k[3] = h * x3p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
-
-        x1[p + 1] = x1[p] + (1 / 6) * (i[0] + (2 * i[1]) + (2 * i[2]) + i[3])
-        x2[p + 1] = x2[p] + (1 / 6) * (j[0] + (2 * j[1]) + (2 * j[2]) + j[3])
-        x3[p + 1] = x3[p] + (1 / 6) * (k[0] + (2 * k[1]) + (2 * k[2]) + k[3])
-
-        t[p + 1] = t[p] + h
-
-        p += 1  # advance flag variable
-
-    if whichy in (0, 1, 2, 3):
-        if whichx == 1:
-            return x1[whichy]
-        elif whichx == 2:
-            return x2[whichy]
-        elif whichx == 3:
-            return x3[whichy]
-        else:
-            return
-
-    return [x1, x2, x3]
-
-
-def plot(a, b, n):  # lower bound, upper bound, num iterations
-    # res = rk(a, b, n)
-    res = adams(a, b, n)
-
-    # Plotting
-    fig = plt.figure()
-    ax = fig.gca(projection='3d')
-    ax.plot(res[0], res[1], res[2], linewidth=1.0, color="darkblue")
-    ax.set_xlabel("X")
-    ax.set_ylabel("Y")
-    ax.set_zlabel("Z")
-    plt.title("Adams-Bashforth-Moutlon Estimation")
-    plt.style.use('ggplot')
-    plt.show()
+class ABM:
+    def __init__(
+        self,
+        f1,
+        f2,
+        f3,
+        t0: float,
+        x0: float,
+        y0: float,
+        z0: float,
+        a: float,
+        b: float,
+        n: int,
+    ):
+        # function 1, func2, func3, initial t0, x0, y0, z0, lower bound, upper, num iterations
+        self.f1 = f1
+        self.f2 = f2
+        self.f3 = f3
+        self.t0 = t0
+        self.x0 = x0
+        self.y0 = y0
+        self.z0 = z0
+        self.a = a
+        self.b = b
+        self.n = n
+
+    def abm_comp(self):
+        self.adams()
+
+    # derivative functions
+    def x1p(self, x1: float, x2: float, x3: float) -> float:
+        return self.f1(x1, x2, x3)
+
+    def x2p(self, x1: float, x2: float, x3: float) -> float:
+        return self.f2(x1, x2, x3)
+
+    def x3p(self, x1: float, x2: float, x3: float) -> float:
+        return self.f3(x1, x2, x3)
+
+    def adams(self) -> list:  # A-B-M computation, using RK4 to start
+        a = self.a
+        b = self.b
+        n = self.n
+
+        h = (b - a) / n
+        y1 = np.zeros(n)  # initialize y1, y2, y3 matrices (and y_predictors)
+        y2 = np.zeros(n)
+        y3 = np.zeros(n)
+
+        y_pred1 = np.zeros(n)
+        y_pred2 = np.zeros(n)
+        y_pred3 = np.zeros(n)
+
+        for E in (0, 1, 2, 3):
+            y1[E] = self.rk(E, 1)  # initial estimations to start A-B-M
+            y2[E] = self.rk(E, 2)
+            y3[E] = self.rk(E, 3)
+
+        for i in range(
+            4, n
+        ):  # Equations change slightly, hard to refactor to make clean
+
+            # x1p
+            y_pred1[i] = y1[i - 1] + (
+                h
+                / 24
+                * (
+                    55 * self.x1p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 59 * self.x1p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + 37 * self.x1p(y1[i - 3], y2[i - 3], y3[i - 3])
+                    - 9 * 0
+                )
+            )
+
+            # x1p
+            y1[i] = y1[i - 1] + (
+                h
+                / 24
+                * (
+                    9 * self.x1p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1])
+                    + 19 * self.x1p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 5 * self.x1p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + self.x1p(y1[i - 3], y2[i - 3], y3[i - 3])
+                )
+            )
+
+            # x2p
+            y_pred2[i] = y2[i - 1] + (
+                h
+                / 24
+                * (
+                    55 * self.x2p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 59 * self.x2p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + 37 * self.x2p(y1[i - 3], y2[i - 3], y3[i - 3])
+                    - 9 * (15 * (-8) - 15)
+                )
+            )
+
+            # x2p
+            y2[i] = y2[i - 1] + (
+                h
+                / 24
+                * (
+                    9 * self.x2p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1])
+                    + 19 * self.x2p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 5 * self.x2p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + self.x2p(y1[i - 3], y2[i - 3], y3[i - 3])
+                )
+            )
+
+            # x3p
+            y_pred3[i] = y3[i - 1] + (
+                h
+                / 24
+                * (
+                    55 * self.x3p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 59 * self.x3p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + 37 * self.x3p(y1[i - 3], y2[i - 3], y3[i - 3])
+                    - 9 * ((15 * 15) - ((8 / 3) * 36))
+                )
+            )
+
+            # x3p
+            y3[i] = y3[i - 1] + (
+                h
+                / 24
+                * (
+                    9 * self.x3p(y_pred1[i - 1], y_pred2[i - 1], y_pred3[i - 1])
+                    + 19 * self.x3p(y1[i - 1], y2[i - 1], y3[i - 1])
+                    - 5 * self.x3p(y1[i - 2], y2[i - 2], y3[i - 2])
+                    + self.x3p(y1[i - 3], y2[i - 3], y3[i - 3])
+                )
+            )
+
+        return [y1, y2, y3]
+
+    def rk(
+        self, which_y: int = 9, which_x: int = 9
+    ) -> list:  # Runge-Kutta 4 computation by itself
+        # which_y indicator for which Y1, Y2, Y3 val is needed to be return by function
+        # which_x means x1, x2, x3 when calling for values Y1, Y2, Y3, etc...
+
+        a = self.a
+        b = self.b
+        n = self.n
+
+        h = (b - a) / n  # step size
+        p = 0  # flag for while loop
+
+        # initialize lists for plot
+        t = np.zeros(int(n + 1))
+        x1 = np.zeros(int(n + 1))
+        x2 = np.zeros(int(n + 1))
+        x3 = np.zeros(int(n + 1))
+
+        # define initial values
+        t[0] = self.t0
+        x1[0] = self.x0
+        x2[0] = self.y0
+        x3[0] = self.z0
+
+        # initialize slots for values of RK4
+        i = np.zeros(4)  # for x1
+        j = np.zeros(4)  # for x2
+        k = np.zeros(4)  # for x3
+
+        # iterate n times
+        while p < n:
+            i[0] = h * self.x1p(x1[p], x2[p], x3[p])
+            j[0] = h * self.x2p(x1[p], x2[p], x3[p])
+            k[0] = h * self.x3p(x1[p], x2[p], x3[p])
+
+            i[1] = h * self.x1p(
+                x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0]
+            )
+            j[1] = h * self.x2p(
+                x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0]
+            )
+            k[1] = h * self.x3p(
+                x1[p] + (1 / 2) * i[0], x2[p] + (1 / 2) * j[0], x3[p] + (1 / 2) * k[0]
+            )
+
+            i[2] = h * self.x1p(
+                x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1]
+            )
+            j[2] = h * self.x2p(
+                x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1]
+            )
+            k[2] = h * self.x3p(
+                x1[p] + (1 / 2) * i[1], x2[p] + (1 / 2) * j[1], x3[p] + (1 / 2) * k[1]
+            )
+
+            i[3] = h * self.x1p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
+            j[3] = h * self.x2p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
+            k[3] = h * self.x3p(x1[p] + i[2], x2[p] + j[2], x3[p] + k[2])
+
+            x1[p + 1] = x1[p] + (1 / 6) * (i[0] + (2 * i[1]) + (2 * i[2]) + i[3])
+            x2[p + 1] = x2[p] + (1 / 6) * (j[0] + (2 * j[1]) + (2 * j[2]) + j[3])
+            x3[p + 1] = x3[p] + (1 / 6) * (k[0] + (2 * k[1]) + (2 * k[2]) + k[3])
+
+            t[p + 1] = t[p] + h
+
+            p += 1  # advance flag variable
+
+        if which_y in (0, 1, 2, 3):
+            if which_x == 1:
+                return x1[which_y]
+            elif which_x == 2:
+                return x2[which_y]
+            elif which_x == 3:
+                return x3[which_y]
+            else:
+                return []
+
+        return [x1, x2, x3]
+
+    def plot3d(self):
+        # res = self.rk()  # to compare to RK solution only
+        res = self.adams()
+
+        # Plotting
+        fig = plt.figure()
+        ax = fig.gca(projection="3d")
+        ax.plot(res[0], res[1], res[2], linewidth=1.0, color="darkblue")
+        ax.set_xlabel("X")
+        ax.set_ylabel("Y")
+        ax.set_zlabel("Z")
+        plt.title("Adams-Bashforth-Moulton Estimation")
+        plt.style.use("ggplot")
+        plt.show()
diff --git a/computepac/euler.py b/computepac/euler.py
index b9acc12..1531fb3 100644
--- a/computepac/euler.py
+++ b/computepac/euler.py
@@ -22,13 +22,15 @@
 import time
 
 
-def euler_forward(f, x, y, a, b, n) -> list:    # function, x0, y0, lower bound, upper bound, num steps
-    x1 = sym.Symbol('x1')
-    x2 = sym.Symbol('x2')
+def euler_forward(
+    f, x, y, a, b, n
+) -> list:  # function, x0, y0, lower bound, upper bound, num steps
+    x1 = sym.Symbol("x1")
+    x2 = sym.Symbol("x2")
     start_time = time.time()
 
-    p = 0                               # flag
-    h = (b - a) / n                     # step size
+    p = 0  # flag
+    h = (b - a) / n  # step size
 
     x_val = []
     y_val = []
@@ -48,13 +50,15 @@ def euler_forward(f, x, y, a, b, n) -> list:    # function, x0, y0, lower bound,
     return [y_val, total_time]
 
 
-def euler_backward(f, x, y, a, b, n) -> list:   # function, x0, y0, lower bound, upper bound, num steps
-    x1 = sym.Symbol('x1')
-    x2 = sym.Symbol('x2')
+def euler_backward(
+    f, x, y, a, b, n
+) -> list:  # function, x0, y0, lower bound, upper bound, num steps
+    x1 = sym.Symbol("x1")
+    x2 = sym.Symbol("x2")
     start_time = time.time()
 
-    p = 0                               # flag
-    h = (b - a) / n                     # step size
+    p = 0  # flag
+    h = (b - a) / n  # step size
 
     x_val = []
     y_val = []
@@ -65,7 +69,7 @@ def euler_backward(f, x, y, a, b, n) -> list:   # function, x0, y0, lower bound,
 
         x += h
         y_p = y + h * f.evalf(subs={x1: x, x2: y})
-        y = y + h * f.evalf(subs={x1: x+h, x2: y_p})
+        y = y + h * f.evalf(subs={x1: x + h, x2: y_p})
 
         p += 1
 
diff --git a/computepac/newton_raphson.py b/computepac/newton_raphson.py
index 210cb27..e010d1c 100644
--- a/computepac/newton_raphson.py
+++ b/computepac/newton_raphson.py
@@ -18,13 +18,13 @@
 
 
 def newton_raphson(f, g, acc) -> list:
-    x = sym.Symbol('x')             # set up x as symbol for functions below
+    x = sym.Symbol("x")  # set up x as symbol for functions below
     start_time = time.time()
-    f_prime = sym.diff(f)           # derivative of inputted function
-    err = []                        # to store errors
-    made_accuracy = 0               # (flag) program has not achieved proper accuracy yet
+    f_prime = sym.diff(f)  # derivative of inputted function
+    err = []  # to store errors
+    made_accuracy = 0  # (flag) program has not achieved proper accuracy yet
 
-    while made_accuracy == 0:       # if desired accuracy not reached yet
+    while made_accuracy == 0:  # if desired accuracy not reached yet
 
         g_1 = g - (f.evalf(subs={x: g}) / f_prime.evalf(subs={x: g}))
 
@@ -32,12 +32,14 @@ def newton_raphson(f, g, acc) -> list:
         err.append(error)
 
         if error < acc:
-            made_accuracy = 1       # signal flag to end computation if accuracy has been reached
+            made_accuracy = (
+                1  # signal flag to end computation if accuracy has been reached
+            )
         else:
-            g = g_1                 # otherwise, continue iterating
+            g = g_1  # otherwise, continue iterating
 
     end_time = time.time()
-    total_time = end_time-start_time
+    total_time = end_time - start_time
     num_itr = len(err)
     root = g
 
diff --git a/computepac/runge_kutta_four.py b/computepac/runge_kutta_four.py
index 5a82342..7a60744 100644
--- a/computepac/runge_kutta_four.py
+++ b/computepac/runge_kutta_four.py
@@ -18,13 +18,15 @@
 import time
 
 
-def rkfour(f, x, y, a, b, n) -> list:  # function, x0, y0, lower bound, upper bound, num steps
-    x1 = sym.Symbol('x1')
-    x2 = sym.Symbol('x2')
-    start_time = time.time() 
+def rkfour(
+    f, x, y, a, b, n
+) -> list:  # function, x0, y0, lower bound, upper bound, num steps
+    x1 = sym.Symbol("x1")
+    x2 = sym.Symbol("x2")
+    start_time = time.time()
 
-    p = 0                       # flag
-    h = (b - a) / n             # step size
+    p = 0  # flag
+    h = (b - a) / n  # step size
 
     x_val = []
     y_val = []
@@ -43,8 +45,8 @@ def rkfour(f, x, y, a, b, n) -> list:  # function, x0, y0, lower bound, upper bo
         y = y + ((1 / 6) * (k0 + (2 * k1) + (2 * k2) + k3))
 
         p += 1
-        
+
     end_time = time.time()
-    total_time = end_time-start_time
+    total_time = end_time - start_time
 
     return [y_val, total_time]
diff --git a/computepac/runge_kutta_two.py b/computepac/runge_kutta_two.py
index 6512cee..74c02df 100644
--- a/computepac/runge_kutta_two.py
+++ b/computepac/runge_kutta_two.py
@@ -19,13 +19,15 @@
 
 
 # Note this is using Heun's constants
-def rktwo(f, x, y, a, b, n) -> list:    # function, x0, y0, lower bound, upper bound, num steps
-    x1 = sym.Symbol('x1')
-    x2 = sym.Symbol('x2')
+def rktwo(
+    f, x, y, a, b, n
+) -> list:  # function, x0, y0, lower bound, upper bound, num steps
+    x1 = sym.Symbol("x1")
+    x2 = sym.Symbol("x2")
     start_time = time.time()
 
-    p = 0                       # flag
-    h = (b - a) / n             # step size
+    p = 0  # flag
+    h = (b - a) / n  # step size
 
     x_val = []
     y_val = []
diff --git a/computepac/utils/utilities.py b/computepac/utils/utilities.py
index 45b4739..1e633f0 100644
--- a/computepac/utils/utilities.py
+++ b/computepac/utils/utilities.py
@@ -4,7 +4,7 @@ def __init__(self):
 
     @staticmethod
     def calc_relative_error(calculated_val: float, correct_val: float) -> float:
-        error = ((correct_val - calculated_val)/correct_val)*100
+        error = ((correct_val - calculated_val) / correct_val) * 100
         return error
 
     @staticmethod
diff --git a/examples/example_convergence.py b/examples/example_convergence.py
index de0296b..c06a2d9 100644
--- a/examples/example_convergence.py
+++ b/examples/example_convergence.py
@@ -36,9 +36,9 @@ def convergence(f, acc):
     plot(f)
 
 
-if __name__ == '__main__':
-    x = Symbol('x')
-    f = Function('f')
+if __name__ == "__main__":
+    x = Symbol("x")
+    f = Function("f")
 
     f = 6 * x ** 3 + 4 * x ** 2 + x + 1  # Enter function to find roots of
 
diff --git a/examples/example_dynamic_Lorenz.py b/examples/example_dynamic_Lorenz.py
index 3c9bddf..862bc8f 100644
--- a/examples/example_dynamic_Lorenz.py
+++ b/examples/example_dynamic_Lorenz.py
@@ -1,9 +1,8 @@
 # Example of the Lorenz Attractor solved by Adams-Bashforth-Moulton
 # in the dynamic_triple module
 #
-# Should be run in an IDE like Anaconda.  May not work with PyCharm.
 
-from computepac import *
+from computepac.dynamic_triple import ABM
 
 
 # Lorenz Attractor equations
@@ -16,8 +15,8 @@ def f2(x, y, z):
 
 
 def f3(x, y, z):
-    return (x * y) - (8/3)*z
+    return (x * y) - (8 / 3) * z
 
 
-abm_comp(f1, f2, f3, 0, 15, 15, 36, 0, 20, 1000)
-
+res = ABM(f1, f2, f3, 0, 15, 15, 36, 0, 20, 1000)
+res.plot3d()
diff --git a/examples/example_euler.py b/examples/example_euler.py
index 259e088..dd53d25 100644
--- a/examples/example_euler.py
+++ b/examples/example_euler.py
@@ -22,22 +22,40 @@
 
 def example_euler(n):
 
-    f = sym.Function('f')
-    x1 = sym.Symbol('x1')
-    x2 = sym.Symbol('x2')
+    f = sym.Function("f")
+    x1 = sym.Symbol("x1")
+    x2 = sym.Symbol("x2")
     f = 6 - 2 * (x2 / x1)
 
     resultf = euler_forward(f, 3, 1, 3, 12, n)
     resultb = euler_backward(f, 3, 1, 3, 12, n)
 
-    print("\nResult of", n, "iterations FORWARD method to estimate f(12) =", resultf[0][n-1])
-    print("Result of", n, "iterations BACKWARD method to estimate f(12) =", resultb[0][n - 1])
+    print(
+        "\nResult of",
+        n,
+        "iterations FORWARD method to estimate f(12) =",
+        resultf[0][n - 1],
+    )
+    print(
+        "Result of",
+        n,
+        "iterations BACKWARD method to estimate f(12) =",
+        resultb[0][n - 1],
+    )
     print("\nTime for FORWARD computation", round(resultf[1], 3), "seconds")
     print("Time for BACKWARD computation", round(resultb[1], 3), "seconds")
     print("\nExact solution known as 23.6875")
-    print("\nFORWARD Relative error:", round(((abs(23.6875 - resultf[0][n-1]) / 23.6875)*100), 5), "%")
-    print("BACKWARD Relative error:", round(((abs(23.6875 - resultb[0][n - 1]) / 23.6875) * 100), 5), "%\n")
-
-
-if __name__ == '__main__':
+    print(
+        "\nFORWARD Relative error:",
+        round(((abs(23.6875 - resultf[0][n - 1]) / 23.6875) * 100), 5),
+        "%",
+    )
+    print(
+        "BACKWARD Relative error:",
+        round(((abs(23.6875 - resultb[0][n - 1]) / 23.6875) * 100), 5),
+        "%\n",
+    )
+
+
+if __name__ == "__main__":
     example_euler(10000)
diff --git a/requirements.txt b/requirements.txt
index 4c7a684..7005f20 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,17 @@
-matplotlib>=3.4.1
-numpy>=1.20.2
-sympy>=1.7.1
-scipy>=1.6.2
+cycler==0.10.0
+iniconfig==1.1.1
+kiwisolver==1.3.1
+matplotlib==3.4.1
+mpmath==1.2.1
+numpy==1.20.2
+packaging==20.9
+Pillow==8.2.0
+pluggy==0.13.1
+py==1.10.0
+pyparsing==2.4.7
+pytest==6.2.3
+python-dateutil==2.8.1
+scipy==1.6.2
+six==1.15.0
+sympy==1.8
+toml==0.10.2
diff --git a/setup.py b/setup.py
index 5c90d7d..34cbdf0 100644
--- a/setup.py
+++ b/setup.py
@@ -3,23 +3,15 @@
 from setuptools import setup, find_packages
 
 setup(
-    name='computepac',
-    version='1.0',
-    description='A Python package of mathematical functions, iterative solvers, numerical analysis, and more.',
-    url='https://github.com/mathemacode/compute-pac',
-    author='mathemacode',
-    classifiers=[
-        'License :: GNU GPLv3'
-        'Programming Language :: Python 3.7'
-    ],
-    keywords='computational mathematics engineering numerical analysis finite difference iterative schemes',
-    package_dir={'': 'computepac'},
-    packages=find_packages(where='computepac'),
-    python_requires='>3.5',
-    install_requires=[
-        'numpy',
-        'sympy',
-        'pytest'
-    ],
-
+    name="computepac",
+    version="1.0",
+    description="A Python package of mathematical functions, iterative solvers, numerical analysis, and more.",
+    url="https://github.com/mathemacode/compute-pac",
+    author="mathemacode",
+    classifiers=["License :: GNU GPLv3" "Programming Language :: Python 3.7"],
+    keywords="computational mathematics engineering numerical analysis finite difference iterative schemes",
+    package_dir={"": "computepac"},
+    packages=find_packages(where="computepac"),
+    python_requires=">3.5",
+    install_requires=["numpy", "sympy", "pytest"],
 )
diff --git a/tests/test_abm.py b/tests/test_abm.py
new file mode 100644
index 0000000..26ae2c1
--- /dev/null
+++ b/tests/test_abm.py
@@ -0,0 +1,19 @@
+from computepac.utils.utilities import Utilities
+from computepac.dynamic_triple import ABM
+
+
+class Test:
+    def test_abm(self):
+        # Lorenz Attractor equations
+        def f1(x, y, z):
+            return 10 * (y - x)
+
+        def f2(x, y, z):
+            return x * (28 - z) - y
+
+        def f3(x, y, z):
+            return (x * y) - (8 / 3) * z
+
+        compute_test = ABM(f1, f2, f3, 0, 15, 15, 36, 0, 20, 1000)
+        ans = compute_test.adams()[2][100]
+        assert Utilities.calc_relative_error(ans, 29.386854067862377) < 0.000001
diff --git a/tests/test_bisection.py b/tests/test_bisection.py
index 01b21f3..730623c 100644
--- a/tests/test_bisection.py
+++ b/tests/test_bisection.py
@@ -5,8 +5,8 @@
 
 class Test:
     def test_bisection(self):
-        f = sym.Function('f')
-        x = sym.Symbol('x')
+        f = sym.Function("f")
+        x = sym.Symbol("x")
         f = (x ** 3) + 2 * (x ** 2) + 10 * x - 20
         ans = bisection(f, 1.0, 2.0, 0.5e-12)[0]
         assert Utilities.check_same_float(ans, 41.0, 1) is True
diff --git a/tests/test_euler.py b/tests/test_euler.py
index 77adfe1..c5d0167 100644
--- a/tests/test_euler.py
+++ b/tests/test_euler.py
@@ -6,18 +6,17 @@
 
 class Test:
     def test_euler_forward(self):
-        f = sym.Function('f')
-        x1 = sym.Symbol('x1')
-        x2 = sym.Symbol('x2')
-        f = 6 - 2*(x2 / x1)
+        f = sym.Function("f")
+        x1 = sym.Symbol("x1")
+        x2 = sym.Symbol("x2")
+        f = 6 - 2 * (x2 / x1)
         ans = euler_forward(f, 3, 1, 3, 6, 100)[0][99]
         assert Utilities.calc_relative_error(ans, 10.72684) < 0.0001
 
-
     def test_euler_backward(self):
-        f = sym.Function('f')
-        x1 = sym.Symbol('x1')
-        x2 = sym.Symbol('x2')
+        f = sym.Function("f")
+        x1 = sym.Symbol("x1")
+        x2 = sym.Symbol("x2")
         f = 6 - 2 * (x2 / x1)
         ans = euler_backward(f, 3, 1, 3, 6, 100)[0][99]
         assert Utilities.calc_relative_error(ans, 10.69012) < 0.0001
diff --git a/tests/test_newton_raphson.py b/tests/test_newton_raphson.py
index f3a8aef..4701f47 100644
--- a/tests/test_newton_raphson.py
+++ b/tests/test_newton_raphson.py
@@ -5,8 +5,8 @@
 
 class Test:
     def test_newton_raphson(self):
-        f = sym.Function('f')
-        x = sym.Symbol('x')
+        f = sym.Function("f")
+        x = sym.Symbol("x")
         f = (x ** 3) + 2 * (x ** 2) + 10 * x - 20
         ans = newton_raphson(f, 2, 0.5e-12)[0]
         assert Utilities.check_same_float(ans, 6.0, 1) is True
diff --git a/tests/test_rkfour.py b/tests/test_rkfour.py
index c7b8cb9..61a68bc 100644
--- a/tests/test_rkfour.py
+++ b/tests/test_rkfour.py
@@ -5,9 +5,9 @@
 
 class Test:
     def test_rkfour(self):
-        f = sym.Function('f')
-        x1 = sym.Symbol('x1')
-        x2 = sym.Symbol('x2')
-        f = (((5 * x1**2) - x2) / sym.exp(x1+x2))
+        f = sym.Function("f")
+        x1 = sym.Symbol("x1")
+        x2 = sym.Symbol("x2")
+        f = ((5 * x1 ** 2) - x2) / sym.exp(x1 + x2)
         ans = round(rkfour(f, 0, 1, 0, 1, 100)[0][99], 5)
         assert Utilities.calc_relative_error(ans, 1.06665) < 0.0001
diff --git a/tests/test_rktwo.py b/tests/test_rktwo.py
index 2e5815e..75ef390 100644
--- a/tests/test_rktwo.py
+++ b/tests/test_rktwo.py
@@ -5,9 +5,9 @@
 
 class Test:
     def test_rktwo(self):
-        f = sym.Function('f')
-        x1 = sym.Symbol('x1')
-        x2 = sym.Symbol('x2')
-        f = (((5 * x1**2) - x2) / sym.exp(x1+x2))
+        f = sym.Function("f")
+        x1 = sym.Symbol("x1")
+        x2 = sym.Symbol("x2")
+        f = ((5 * x1 ** 2) - x2) / sym.exp(x1 + x2)
         ans = round(rktwo(f, 0, 1, 0, 1, 100)[0][99], 5)
         assert Utilities.calc_relative_error(ans, 1.00058) < 0.0001
diff --git a/tests/test_sample.py b/tests/test_sample.py
index be2265d..9cda1b4 100644
--- a/tests/test_sample.py
+++ b/tests/test_sample.py
@@ -2,7 +2,7 @@
 
 
 def mult_zero(x):
-    return x*0
+    return x * 0
 
 
 class Test: