diff --git a/CITATION.cff b/CITATION.cff
index e77905a..f4df2c1 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -1,4 +1,4 @@
-cff-version: 1.1.0
+cff-version: 1.2.0
 message: "If you use this software, please cite it as below."
 title: Gaspype
 abstract: Gaspype is a performant library for thermodynamic calculations with ideal gases
@@ -9,11 +9,11 @@ authors:
     affiliation: "German Aerospace Center (DLR)"
     address: "Linder Höhe"
     city: Köln
-version: v1.1.3
+version: v1.1.4
 date-released: "2025-06-24"
-#identifiers:
-#  - description: This is the collection of archived snapshots of all versions of Gaspype
-#    type: doi
-#    value: ""
+identifiers:
+  - description: This is the collection of archived snapshots of all versions of Gaspype
+    type: doi
+    value: "10.5281/zenodo.17047601"
 license: MIT
 repository-code: "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"
diff --git a/examples/methane_mixtures.md b/examples/methane_mixtures.md
index 2fa8448..2955892 100644
--- a/examples/methane_mixtures.md
+++ b/examples/methane_mixtures.md
@@ -37,7 +37,7 @@ ax.set_ylabel("molar fraction")
 ax.set_ylim(0, 1.1)
 #ax.set_xlim(0, 100)
 ax.plot(ratio, equilibrium_h2o.get_x())
-ax.legend(fs.active_species)
+ax.legend(fs.species)
 ```
 
 Equilibrium calculation for methane CO2 mixtures:
@@ -56,5 +56,5 @@ ax.set_ylabel("molar fraction")
 ax.set_ylim(0, 1.1)
 #ax.set_xlim(0, 100)
 ax.plot(ratio, equilibrium_co2.get_x())
-ax.legend(fs.active_species)
+ax.legend(fs.species)
 ```
diff --git a/pyproject.toml b/pyproject.toml
index 8ba3700..6424fa1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "gaspype"
-version = "1.1.3"
+version = "1.1.4"
 authors = [
   { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
 ]
@@ -14,7 +14,6 @@ classifiers = [
 ]
 dependencies = [
   "numpy>2.0.0",
-  "scipy>1.12.0",
 ]
 
 [project.urls]
@@ -42,7 +41,6 @@ dev = [
     "cantera",
     "pyyaml>=6.0.1",
     "types-PyYAML",
-    "scipy-stubs",
     "matplotlib"
 ]
 doc_build = [
diff --git a/src/gaspype/_main.py b/src/gaspype/_main.py
index 1bd0afe..33f96d9 100644
--- a/src/gaspype/_main.py
+++ b/src/gaspype/_main.py
@@ -2,7 +2,7 @@ import numpy as np
 from numpy.typing import NDArray
 from typing import Sequence, Any, TypeVar, Iterator, overload, Callable
 from math import log as ln, ceil
-from scipy.linalg import null_space
+from ._numerics import null_space
 from gaspype._phys_data import atomic_weights, db_reader
 import re
 import pkgutil
@@ -90,7 +90,7 @@ class fluid_system:
 
         self._t_offset = int(t_min)
         self.species = species
-        self.active_species = species
+        self.active_species = species  # for backward compatibility
         element_compositions: list[dict[str, int]] = list()
 
         for i, s in enumerate(species):
@@ -220,7 +220,7 @@ class fluid:
                 The array can be multidimensional, the size of the last dimension
                 must match the number of species defined for the fluid_system.
                 The indices of the last dimension correspond to the indices in
-                the active_species list of the fluid_system.
+                the species list of the fluid_system.
             fs: Reference to a fluid_system. Is optional if composition is
                 defined by a dict. If not specified a new fluid_system with
                 the components from the dict is created.
@@ -585,7 +585,7 @@ class elements:
                 The array can be multidimensional, the size of the last dimension
                 must match the number of elements used in the fluid_system.
                 The indices of the last dimension correspond to the indices in
-                the active_species list of the fluid_system.
+                the species list of the fluid_system.
             fs: Reference to a fluid_system.
             shape: Tuple or list for the dimensions the fluid array. Can
                 only be used if composition argument is a dict. Otherwise
diff --git a/src/gaspype/_numerics.py b/src/gaspype/_numerics.py
new file mode 100644
index 0000000..dc4717a
--- /dev/null
+++ b/src/gaspype/_numerics.py
@@ -0,0 +1,21 @@
+import numpy as np
+from .typing import FloatArray
+
+
+def null_space(A: FloatArray) -> FloatArray:
+    """
+    Compute an orthonormal basis for the null space of A using NumPy SVD.
+
+    Args:
+        A: Input matrix of shape (m, n)
+
+    Return:
+        Null space vectors as columns, shape (n, n - rank)
+    """
+    u, s, vh = np.linalg.svd(A, full_matrices=True)
+    M, N = u.shape[0], vh.shape[1]
+    rcond = np.finfo(s.dtype).eps * max(M, N)
+    tol = np.amax(s, initial=0.) * rcond
+    num = np.sum(s > tol, dtype=int)
+    Q = vh[num:, :].T.conj()
+    return Q
diff --git a/src/gaspype/_solver.py b/src/gaspype/_solver.py
index d3f7c96..100aa94 100644
--- a/src/gaspype/_solver.py
+++ b/src/gaspype/_solver.py
@@ -1,10 +1,19 @@
-from typing import Literal, Any
-from scipy.optimize import minimize, root
+from typing import Literal, Any, TYPE_CHECKING
 import numpy as np
 from ._main import elements, fluid, fluid_system
 from .typing import NDFloat, FloatArray
 from .constants import p0, epsy
 
+if TYPE_CHECKING:
+    def minimize(*a: Any, **b: Any) -> dict[str, FloatArray]:
+        ...
+else:
+    try:
+        from scipy.optimize import minimize
+    except ImportError:
+        def minimize(*a, **b):
+            raise ImportError('scipy is required for the "gibs minimization" solver')
+
 
 def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
     """
@@ -13,8 +22,7 @@ def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> N
     Solvers:
         - **system of equations** (default): Finds the root for a system of
           equations covering a minimal set of equilibrium equations and elemental balance.
-          The minimal set of equilibrium equations is derived by SVD using the null_space
-          implementation of scipy.
+          The minimal set of equilibrium equations is derived by SVD calculating null space.
 
         - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
           the elemental composition constant using the SLSQP implementation of scipy
@@ -60,7 +68,7 @@ def equilibrium_gmin(fs: fluid_system, element_composition: FloatArray, t: float
 
     start_composition_array = np.ones_like(fs.species, dtype=float)
     sol = np.array(minimize(gibbs_rt, start_composition_array, args=(grt, p_rel), method='SLSQP',
-                   bounds=bnds, constraints=cons, options={'maxiter': 2000, 'ftol': 1e-12})['x'], dtype=NDFloat)  # type: ignore
+                   bounds=bnds, constraints=cons, options={'maxiter': 2000, 'ftol': 1e-12})['x'], dtype=NDFloat)
 
     return sol
 
@@ -72,7 +80,6 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
     element_norm_log = np.log(element_norm + epsy)
 
     a = fs.array_stoichiometric_coefficients
-    a_sum = np.sum(a)
     el_matrix = fs.array_species_elements.T
 
     # Log equilibrium constants for each reaction equation
@@ -82,42 +89,49 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
     bp = b - np.sum(a * np.log(p / p0), axis=1)
 
     # Calculating the maximum possible amount for each species based on the elements
-    species_max = np.min(element_norm / (fs.array_species_elements + epsy), axis=1)
-    logn_start = np.log(species_max + epsy)
+    species_max = np.min((element_norm + epsy) / (fs.array_species_elements + epsy), axis=1)
+    species_max_log = np.log(species_max + epsy)
 
-    # global count
-    # count = 0
-
-    weighting = 100
+    # Prepare constant arrays
+    j_eq_eye = np.eye(len(species_max))
+    j_eq_ones = np.ones((len(species_max), 1))
 
     def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
-        # global count
-        # count += 1
-        # print('------', count)
-        # assert count < 100
-        n = np.exp(logn)
+        n: FloatArray = np.exp(logn)  # n is the molar amount normalized by el_max
         n_sum = np.sum(n)
 
         # Residuals from equilibrium equations:
-        resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
+        resid_eq = a @ (logn - np.log(n_sum)) - bp
 
-        # Jacobian:
-        j_eq = a - a_sum * n / n_sum
+        # Jacobian for equilibrium equations:
+        j_eq = a @ (j_eq_eye - j_eq_ones * n / np.sum(n))
 
         # Residuals from elemental balance:
-        el_sum = np.dot(el_matrix, n)
-        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
+        el_sum_norm = np.dot(el_matrix, n) + epsy
+        resid_ab = np.log(el_sum_norm) - element_norm_log
+        #print(f'* resid_eq: {resid_eq}      resid_ab: {resid_ab}            {element_norm}')
 
-        # print(el_sum, element_norm)
-
-        # Jacobian
-        j_ab = weighting * el_matrix * n / el_sum[:, np.newaxis]
+        # Jacobian for elemental balance:
+        j_ab = el_matrix * n / el_sum_norm[:, None]
 
         return (np.hstack([resid_eq, resid_ab]), np.concatenate([j_eq, j_ab], axis=0))
 
-    ret = root(residuals, logn_start, jac=True, tol=1e-10)
-    n = np.exp(np.array(ret['x'], dtype=NDFloat))
-    # print(ret)
+    logn: FloatArray = species_max_log  # Set start values
+
+    for i in range(30):
+        rF, J = residuals(logn)
+
+        delta = np.linalg.solve(J, -rF)
+
+        logn = logn + delta
+        logn = np.minimum(logn, species_max_log + 1)
+
+        #print(f'{i} F: {np.linalg.norm(rF):.5f} lognmin={np.min(logn):.3f}, lognmax={np.max(logn):.3f}, delta={np.linalg.norm(delta):.3f} logn=')
+        if np.linalg.norm(rF) < 1e-10:
+            #print(f'Converged in {i} iterations')
+            break
+
+    n = np.exp(logn)
 
     return n * el_max
 
diff --git a/tests/benchmark_equalibrium.py b/tests/benchmark_equalibrium.py
index 748e98d..4c0fe3a 100644
--- a/tests/benchmark_equalibrium.py
+++ b/tests/benchmark_equalibrium.py
@@ -44,6 +44,10 @@ eq_gaspype = gp.equilibrium(fluid, t=temperatures, p=pressure)
 elapsed_gaspype = time.perf_counter() - t0
 print(f"Gaspype: {elapsed_gaspype:.4f} s")
 
+# Check if elemental balance of result is correct
+el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2)
+assert np.all(el_err < 1e-20)
+
 # -----------------------
 # Compare first 5 results
 # -----------------------
diff --git a/tests/test_equalibrium.py b/tests/test_equalibrium.py
new file mode 100644
index 0000000..26d9eb3
--- /dev/null
+++ b/tests/test_equalibrium.py
@@ -0,0 +1,25 @@
+import gaspype as gp
+
+
+def test_single_equilibrium():
+    # Compare equilibrium calculations to Cantera results
+
+    # gp.set_solver('system of equations')
+    # gp.set_solver('gibs minimization')
+
+    # fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+    fs = gp.fluid_system(['CH4', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+    # fs = gp.fluid_system([s for s in flow1.species_names if s in gps])
+
+    composition = gp.elements({'H': 2, 'O': 0, 'C': 0}, fs)
+
+    t = 1495 + 273.15  # K
+    p = 1e5  # Pa
+
+    fl = gp.equilibrium(composition, t, p)
+
+    print(fl)
+
+
+if __name__ == "__main__":
+    test_single_equilibrium()
diff --git a/tests/test_results.py b/tests/test_results.py
index 165f21a..50a4f9c 100644
--- a/tests/test_results.py
+++ b/tests/test_results.py
@@ -55,6 +55,7 @@ def test_nh3_data():
 def test_equilibrium():
     # Compare equilibrium calculations to Cycle-Tempo results
     df = pd.read_csv('tests/test_data/cycle_temp_matlab_ref.csv', sep=';', decimal=',').fillna(0)
+    #gp.set_solver('gibs minimization')
     fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'C4H10,n-butane', 'H2O', 'H2', 'CO2', 'CO'])
 
     for index, row in df.iterrows():
@@ -84,7 +85,7 @@ def test_equilibrium():
             result_values = gp.equilibrium(fl, t, p).array_fractions
 
             print(index, gp.get_solver(), '----')
-            print(molar_comp)
+            print('Species:        ' + ''.join(f"{s:14}" for s in fs.species))
             outp(result_values, 'Under test: ')
             outp(reference_values, 'Reference:  ')
 
@@ -123,3 +124,7 @@ def test_carbon():
             assert result_values > 0.9
         else:
             assert result_values < 1.1
+
+
+if __name__ == '__main__':
+    test_equilibrium()
diff --git a/tests/test_results_cantera.py b/tests/test_results_cantera.py
index fb397f2..018c6e6 100644
--- a/tests/test_results_cantera.py
+++ b/tests/test_results_cantera.py
@@ -28,14 +28,14 @@ def test_equilibrium_cantera():
 
     flow1 = ct.Solution('gri30.yaml')  # type: ignore
     ct_results = []
-    comp = (composition.array_fractions[i, j] for i in range(composition.shape[0]) for j in range(composition.shape[1]))
+    comp = [composition.array_fractions[i, j] for i in range(composition.shape[0]) for j in range(composition.shape[1])]
     for c in comp:
         comp_dict = {s: v for s, v in zip(fs.species, c)}
 
         flow1.TP = t, p
         flow1.X = comp_dict
         flow1.equilibrate('TP')  # type: ignore
-        indeces = [i for flsn in fs.active_species for i, sn in enumerate(flow1.species_names) if flsn == sn]  # type: ignore
+        indeces = [i for flsn in fs.species for i, sn in enumerate(flow1.species_names) if flsn == sn]  # type: ignore
         ct_results.append(flow1.X[indeces])  # type: ignore
 
         #if flow1.X[indeces][0] > 0.01:
@@ -45,8 +45,33 @@ def test_equilibrium_cantera():
 
     deviations = np.abs(gp_result_array - ct_result_array)
 
-    for dev, gp_comp_result, ct_comp_result in zip(deviations, gp_result_array, ct_result_array):
-        print(f"Composition: {gp_comp_result} / {ct_comp_result}")
+    for dev, gp_comp_result, ct_comp_result, c in zip(deviations, gp_result_array, ct_result_array, comp):
+        comp_dict = {s: v for s, v in zip(fs.species, c)}
+        print(f"Inp. Composition: {comp_dict}")
+        print(f"Res. Composition: {gp_comp_result}")
+        print(f"Ref. Composition: {ct_comp_result}")
+        print("---")
         assert np.all(dev < 0.04), f"Deviateion: {dev}"
 
     assert np.mean(deviations) < 2e-4
+
+
+def test_cantera():
+
+    t = 1495 + 273.15  # K
+    p = 1e5  # Pa
+
+    flow1 = ct.Solution('gri30.yaml')  # type: ignore
+    flow1.TP = t, p
+    inp_comp = {'CH4': 0.0, 'H2O': 0.0, 'H2': 0.9508196721311476, 'CO2': 0.0, 'CO': 0.0, 'O2': 0.04918032786885246}
+    flow1.X = inp_comp
+    flow1.equilibrate('TP')  # type: ignore
+
+    results: dict[str, float] = {sn: float(flow1.X[i]) for flsn in inp_comp for i, sn in enumerate(flow1.species_names) if flsn == sn}  # type: ignore
+
+    print(inp_comp)
+    print(results)
+
+
+if __name__ == "__main__":
+    test_equilibrium_cantera()