diff --git a/Examples/Mesh_tests/Ejemplo_axial.json b/Examples/Mesh_tests/Ejemplo_axial.json
new file mode 100644
index 0000000..e1df8fc
--- /dev/null
+++ b/Examples/Mesh_tests/Ejemplo_axial.json
@@ -0,0 +1 @@
+{"nodes": [[0.0], [1.5], [3.0], [4.0], [5.0]], "dictionary": [[0, 1], [1, 2], [2, 3], [3, 4]], "types": ["L1V", "L1V", "L1V", "L1V"], "regions": [], "ebc": [[0, 0]], "nbc": [[-1, -850]], "nvn": 1, "ngdl": 5, "holes": [], "fillets": [], "solutions": [{"info": {"solver-type": "lineal"}, "U": [[0.0], [-8.116902097686659e-05], [-0.00016233804195373317], [-0.0003126510437627454], [-0.0004629640455717576]]}], "properties": {"verbose": false, "name": "", "problem": "EDO1D", "a": null, "b": null, "c": null, "WARNING": "It's not possible lo save callables", "duration": 0.756277, "description": "FEM problem using the Generic 1D second order diferential equation formulation.,DOF: 5,Elements: 4,Solver: lineal,EBC: 1,NBC: 1"}}
\ No newline at end of file
diff --git a/Examples/example_axial.py b/Examples/example_axial.py
new file mode 100644
index 0000000..996c4df
--- /dev/null
+++ b/Examples/example_axial.py
@@ -0,0 +1,34 @@
+if __name__ == '__main__':
+    from FEM import EDO1D
+    from FEM.Geometry import Geometry1D
+    import matplotlib.pyplot as plt
+    import numpy as np
+    E = 20000000
+    r1 = 0.5
+    r2 = 0.3
+    A1 = r1**2*np.pi
+    A2 = r2**2*np.pi
+
+    def EA(x):
+        EE = (x <= 3)*E*A1
+        EE += (x > 3)*E*A2
+        return EE
+
+    geo = Geometry1D([[0, 1],
+                      [1, 2],
+                      [2, 3],
+                      [3, 4]],
+
+                     [[0.0],
+                      [1.5],
+                      [3],
+                      [4],
+                      [5]],
+
+                     ["L1V"]*4, 1)
+    geo.setCbe([[0, 0]])
+    geo.cbn = [[-1, -850]]
+    O = EDO1D(geo, EA, lambda x: 0, lambda x: 0)
+    O.solve(plot=True)
+    O.exportJSON("Ejemplo_axial.json")
+    plt.show()
diff --git a/TODO.md b/TODO.md
index bf38f48..f68b3a8 100644
--- a/TODO.md
+++ b/TODO.md
@@ -2,7 +2,6 @@
     -   [x] Brick orden 2.
     -   [x] Pirámide orden 1.
     -   [x] Pirámide orden 2.
--   [ ] Los elementos deben agruparse en un solo paquete.
 -   [x] Modificar geometría para que sea una clase mas general.
     -   [x] La clase gemoetría 2D debe tener los métodos que se tienen en la clase geometría actualmente.
     -   [x] Se deben agregar métodos equivalentes en la clase geometría 1D y 3D.
diff --git a/a.py b/a.py
deleted file mode 100644
index 206d430..0000000
--- a/a.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import numpy as np
-import FEM
-from tqdm import tqdm
-
-F1 = FEM.importJSON("Examples/Mesh_tests/AYO_TAPON.json", {}, fast=True)
-F2 = FEM.importJSON("Examples/Mesh_tests/AYO_NO_TAPON.json", {}, fast=True)
-
-
-def error(LDUM, HDUM, nvn=3):
-    errors = []
-    for i in tqdm(range(len(HDUM.geometry.gdls))):
-        coord = HDUM.geometry.gdls[i]
-        potential = LDUM.geometry.KDTree.query(coord, 40, workers=-1)[1]
-        flag = True
-        errores = []
-        for eidx in potential:
-            e = LDUM.elements[eidx]
-            z = e.inverseMapping(coord.reshape([3, 1]))
-            zz = np.sum(z)
-            errores.append(f"{z.flatten()}, {zz}")
-            tol = 1e-6
-            if zz <= (1+tol) and np.all(z.flatten() >= -tol):
-                uLDUMx, uLDUMy, uLDUMz = e.giveSolutionPoint(z)[1]
-                uHDUMx = HDUM.U[nvn*i]
-                uHDUMy = HDUM.U[nvn*i + 1]
-                uHDUMz = HDUM.U[nvn*i + 2]
-                errors.append(abs(uLDUMx-uHDUMx))
-                errors.append(abs(uLDUMy-uHDUMy))
-                errors.append(abs(uLDUMz-uHDUMz))
-                flag = False
-                break
-        if flag:
-            print(f"HEY coord {i}, {coord} NOT WORKED! {errores}")
-            break
-    return np.array(errors)
-
-
-U = error(F1, F2)
-np.savetxt("errores.txt", U)
-print()
diff --git a/src/FEM/Elements/Element.py b/src/FEM/Elements/Element.py
index 12b54af..4c0daee 100644
--- a/src/FEM/Elements/Element.py
+++ b/src/FEM/Elements/Element.py
@@ -46,8 +46,8 @@ def __init__(self, coords: np.ndarray, _coords: np.ndarray, gdl: np.ndarray, bor
         self._x, self._p = self.T(self.Z.T)
         self.jacs, self.dpz = self.J(self.Z.T)
         self._xcenter = self.T(self.center.T)[0].flatten()
-        if not self.border:
-            if not self.fast:
+        if not self.border:  # Border elements do not satisfy coordinate transformation requirements
+            if not self.fast:  # Fast elements are for sparse matrix assembly
                 self.Ke = np.zeros([self.n, self.n])
                 self.Fe = np.zeros([self.n, 1])
                 self.Qe = np.zeros([self.n, 1])
@@ -59,7 +59,7 @@ def __init__(self, coords: np.ndarray, _coords: np.ndarray, gdl: np.ndarray, bor
         self.Ue = np.zeros(self.gdl.shape)
 
     def restartMatrix(self) -> None:
-        """Sets all element matrices and vectors to 0 state
+        """Sets all element matrices and vectors to 0
         """
         if not self.border:
             self.Ke[:, :] = 0.0
diff --git a/src/FEM/Geometry/Region.py b/src/FEM/Geometry/Region.py
index ae080d0..e5f5e0f 100644
--- a/src/FEM/Geometry/Region.py
+++ b/src/FEM/Geometry/Region.py
@@ -76,7 +76,7 @@ class Region2D(Region):
     """
 
     def __init__(self, coords: np.ndarray, **kargs) -> None:
-        """Creates a 2D region using a 2D Element. The number of coordinates definesthe type of element.
+        """Creates a 2D region using a 2D Element. The number of coordinates defines the type of element.
 
         Args:
             coords (np.ndarray): Coordinate matrix. Must be of four rows and 3 columns.