From aa44f69c7c0d375a9cefa8369fd63672025e363a Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 22 Feb 2026 11:08:32 +0100
Subject: [PATCH 01/26] Added capability to find other fields in LHSs and
 layers

---
 layercake/arithmetic/equation.py |  3 +++
 layercake/bakery/layers.py       | 12 ++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 1229991..c405ce2 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -91,6 +91,9 @@ def other_fields(self):
             for term in equation_term.terms:
                 if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                     other_fields.append(term.field)
+        for term in self.lhs_term.terms:
+            if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
+                other_fields.append(term.field)
         return other_fields
 
     @property
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 3ea036d..02a07ea 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -93,6 +93,18 @@ def fields(self):
             fields_list.append(eq.field)
         return fields_list
 
+    @property
+    def other_fields(self):
+        """list(~field.Field): Returns the list of dynamical fields of other layers, i.e. the fields whose time
+        evolution is provided by the partial differential equations of other layers."""
+        layer_fields = self.fields
+        other_fields = list()
+        for eq in self.equations:
+            for field in eq.other_fields:
+                if field not in layer_fields:
+                    other_fields.append(field)
+        return other_fields
+
     @property
     def parameters(self):
         """list(~parameter.Parameter): Returns the list of parameters of the layer, i.e. the explicit parameters

From 75c8a54e20812c084f22a8a52d35df94adda4885 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 22 Feb 2026 11:12:51 +0100
Subject: [PATCH 02/26] Linting

---
 layercake/arithmetic/equation.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index c405ce2..c433cda 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -33,7 +33,7 @@ class Equation(object):
         Must be a single term, possibly a combination
         through :class:`~layercake.arithmetic.terms.base.OperationOnTerms` operations.
     name: str, optional
-        Optional name for the equation.
+        Name for the equation.
 
     Attributes
     ----------
@@ -246,10 +246,10 @@ def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_ch
         else:
             latex_string = r'\frac{\partial}{\partial t} ' + lhs
 
-        fterm = self.terms[0].latex
+        first_term = self.terms[0].latex
         if drop_first_rhs_char:
-            fterm = fterm[2:]
-        latex_string += ' = ' + fterm
+            first_term = first_term[2:]
+        latex_string += ' = ' + first_term
 
         for term in self.terms[1:]:
             latex_string += term.latex

From 67e7655e734bc0c573d01406364c822996d934cf Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 22 Feb 2026 11:42:35 +0100
Subject: [PATCH 03/26] Removing the setting of field in lhs

Seems not needed and harmful if other fields involved
---
 layercake/arithmetic/equation.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index c433cda..00daec6 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -28,7 +28,7 @@ class Equation(object):
     ----------
     field: ~field.Field
         The spatial field over which the partial differential equation.
-    lhs_term: ~arithmetic.terms.base.ArithmeticTerms
+    lhs_term: ~arithmetic.terms.base.ArithmeticTerms or list(~arithmetic.terms.base.ArithmeticTerms)
         Term on the left-hand side of the equation.
         Must be a single term, possibly a combination
         through :class:`~layercake.arithmetic.terms.base.OperationOnTerms` operations.
@@ -55,7 +55,6 @@ def __init__(self, field, lhs_term, name=''):
         self.field._equation = self
         self.terms = list()
         self.lhs_term = lhs_term
-        self.lhs_term.field = self.field
         self.name = name
         self._layer = None
         self._cake = None

From fd39ac4c161aad6b74b01c41527000493d253b07 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sat, 28 Feb 2026 18:34:16 +0100
Subject: [PATCH 04/26] Added the other fields in lhs property for equations

---
 layercake/arithmetic/equation.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 00daec6..625a5cd 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -90,6 +90,12 @@ def other_fields(self):
             for term in equation_term.terms:
                 if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                     other_fields.append(term.field)
+        other_fields = other_fields + self.other_fields_in_lhs
+        return other_fields
+
+    @property
+    def other_fields_in_lhs(self):
+        other_fields = list()
         for term in self.lhs_term.terms:
             if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                 other_fields.append(term.field)

From 928f508d3a86c7db3a353633d06ce0a7af6f81d8 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sat, 28 Feb 2026 19:29:08 +0100
Subject: [PATCH 05/26] Increasing Equation class clarity

and preparing for more complicate lhs
---
 layercake/arithmetic/equation.py | 168 +++++++++++++++++++++++--------
 layercake/bakery/layers.py       |  18 ++--
 2 files changed, 137 insertions(+), 49 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 625a5cd..128563f 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -16,6 +16,8 @@
 from layercake.utils import isin
 
 
+# TODO: deal with the cases where lists of terms are empty - maybe nothing is needed but to check
+
 class Equation(object):
     """Main class to define and specify partial differential equations.
 
@@ -28,10 +30,8 @@ class Equation(object):
     ----------
     field: ~field.Field
         The spatial field over which the partial differential equation.
-    lhs_term: ~arithmetic.terms.base.ArithmeticTerms or list(~arithmetic.terms.base.ArithmeticTerms)
-        Term on the left-hand side of the equation.
-        Must be a single term, possibly a combination
-        through :class:`~layercake.arithmetic.terms.base.OperationOnTerms` operations.
+    lhs_terms: ~arithmetic.terms.base.ArithmeticTerms or list(~arithmetic.terms.base.ArithmeticTerms), optional
+        Terms on the left-hand side of the equation. At least one must involve the field defined above.
     name: str, optional
         Name for the equation.
 
@@ -39,9 +39,9 @@ class Equation(object):
     ----------
     field: ~field.Field
         The spatial field over which the partial differential equation.
-    terms: list(~arithmetic.terms.base.ArithmeticTerms)
+    rhs_terms: list(~arithmetic.terms.base.ArithmeticTerms)
         List of additive terms in the right-hand side of the equation.
-    lhs_term: ~arithmetic.terms.base.ArithmeticTerms
+    lhs_terms: ListOfArithmeticTerms(~arithmetic.terms.base.ArithmeticTerms)
         Term on the left-hand side of the equation.
     name: str
         Optional name for the equation.
@@ -49,12 +49,17 @@ class Equation(object):
 
     _t = Symbol('t')
 
-    def __init__(self, field, lhs_term, name=''):
+    def __init__(self, field, lhs_terms=None, name=''):
 
         self.field = field
         self.field._equation = self
-        self.terms = list()
-        self.lhs_term = lhs_term
+        self.rhs_terms = ListOfAdditiveArithmeticTerms()
+        if lhs_terms is None:
+            self.lhs_terms = ListOfAdditiveArithmeticTerms()
+        elif isinstance(lhs_terms, list):
+            self.lhs_terms = ListOfAdditiveArithmeticTerms(lhs_terms)
+        else:
+            self.lhs_terms = ListOfAdditiveArithmeticTerms([lhs_terms])
         self.name = name
         self._layer = None
         self._cake = None
@@ -82,12 +87,40 @@ def add_rhs_terms(self, terms):
         for t in terms:
             self.add_rhs_term(t)
 
+    def add_lhs_term(self, term):
+        """Add a term to the left-hand side of the equation.
+
+        Parameters
+        ----------
+        term: ~arithmetic.terms.base.ArithmeticTerms
+            Term to be added to the left-hand side of the equation.
+        """
+        if not issubclass(term.__class__, ArithmeticTerms):
+            raise ValueError('Provided term must be a valid ArithmeticTerm object.')
+        self.lhs_terms.append(term)
+
+    def add_lhs_terms(self, terms):
+        """Add multiple terms to the left-hand side of the equation.
+
+        Parameters
+        ----------
+        terms: list(~arithmetic.terms.base.ArithmeticTerms)
+            Terms to be added to the left-hand side of the equation.
+        """
+        for t in terms:
+            self.add_lhs_term(t)
+
+    @property
+    def terms(self):
+        """Alias for the list of RHS arithmetic terms."""
+        return self.terms
+
     @property
     def other_fields(self):
         """list(~field.Field): List of additional fields present in the equation."""
         other_fields = list()
         for equation_term in self.terms:
-            for term in equation_term.terms:
+            for term in equation_term.rhs_terms:
                 if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                     other_fields.append(term.field)
         other_fields = other_fields + self.other_fields_in_lhs
@@ -96,7 +129,7 @@ def other_fields(self):
     @property
     def other_fields_in_lhs(self):
         other_fields = list()
-        for term in self.lhs_term.terms:
+        for term in self.lhs_terms:
             if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                 other_fields.append(term.field)
         return other_fields
@@ -106,10 +139,10 @@ def parameter_fields(self):
         """list(~field.ParameterField): List of non-dynamical parameter fields present in the equation."""
         parameter_fields = list()
         for equation_term in self.terms:
-            for term in equation_term.terms:
+            for term in equation_term.rhs_terms:
                 if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
                     parameter_fields.append(term.field)
-        for term in self.lhs_term.terms:
+        for term in self.lhs_terms:
             if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
                 parameter_fields.append(term.field)
         return parameter_fields
@@ -118,7 +151,7 @@ def parameter_fields(self):
     def parameters(self):
         """list(~parameter.Parameter): List of parameters present in the equation."""
         parameters_list = list()
-        for term in self.terms + [self.lhs_term]:
+        for term in self.terms + self.lhs_terms:
             params_list = term.parameters
             for param in params_list:
                 if not isin(param, parameters_list):
@@ -139,64 +172,49 @@ def parameters_symbols(self):
     @property
     def symbolic_expression(self):
         """~sympy.core.expr.Expr: Symbolic expression of the equation."""
-        rterm = S.Zero
-        for term in self.terms:
-            rterm += term.symbolic_expression
-        return Eq(self.symbolic_lhs.diff(self._t, evaluate=False), rterm)
+        return Eq(self.symbolic_lhs.diff(self._t, evaluate=False), self.rhs_terms.symbolic_expression)
 
     @property
     def numerical_expression(self):
         """~sympy.core.expr.Expr: Expression of the equation with parameters replaced by their
         configured values."""
-        rterm = S.Zero
-        for term in self.terms:
-            rterm += term.numerical_expression
-        return Eq(self.numerical_lhs.diff(self._t, evaluate=False), rterm)
+        return Eq(self.numerical_lhs.diff(self._t, evaluate=False), self.rhs_terms.symbolic_expression)
 
     @property
     def symbolic_rhs(self):
         """~sympy.core.expr.Expr: Symbolic expression of the right-hand side of the equation."""
-        rterm = S.Zero
-        for term in self.terms:
-            rterm += term.symbolic_expression
-        return rterm
+        return self.rhs_terms.symbolic_expression
 
     @property
     def numerical_rhs(self):
         """~sympy.core.expr.Expr: Expression of the right-hand side of the equation with
         parameters replaced by their configured values."""
-        rterm = S.Zero
-        for term in self.terms:
-            rterm += term.numerical_expression
-        return rterm
+        return self.rhs_terms.numerical_expression
 
     @property
     def symbolic_lhs(self):
         """~sympy.core.expr.Expr: Symbolic expression of the left-hand side of the equation."""
-        return self.lhs_term.symbolic_expression
+        return self.lhs_terms.symbolic_expression
 
     @property
     def numerical_lhs(self):
         """~sympy.core.expr.Expr: Expression of the left-hand side of the equation with
         parameters replaced by their configured values."""
-        return self.lhs_term.numerical_expression
+        return self.lhs_terms.numerical_expression
 
     @property
     def lhs_inner_products(self):
-        """~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.tensor.array.ImmutableSparseNDimArray or sparse.COO(float): Left-hand
+        """list(~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.tensor.array.ImmutableSparseNDimArray) or list(sparse.COO(float)): Left-hand
         side inner products of the equation, if available."""
-        return self.lhs_term.inner_products
+        return [term.inner_products for term in self.lhs_terms]
 
     @property
     def maximum_rank(self):
         """int: Maximum rank of the right-hand side terms tensors."""
-        rhs_max_rank = 0
-        for term in self.terms:
-            rhs_max_rank = max(rhs_max_rank, term.rank)
-        return max(rhs_max_rank, self.lhs_term.rank)
+        return max(self.terms.maximum_rank, self.lhs_terms.maximum_rank)
 
     def compute_lhs_inner_products(self, basis, numerical=False, timeout=None, num_threads=None, permute=False):
-        """Compute the inner products tensor of the left-hand side term.
+        """Compute the inner products tensor of the left-hand side terms.
 
         Parameters
         ----------
@@ -219,7 +237,7 @@ def compute_lhs_inner_products(self, basis, numerical=False, timeout=None, num_t
             from 1 to the rank of the tensor.
             Default to `False`, i.e. no permutation is applied.
         """
-        self.lhs_term.compute_inner_products(basis, numerical, timeout, num_threads, permute)
+        self.lhs_terms.compute_inner_products(basis, numerical, timeout, num_threads, permute)
 
     def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_char=False):
         """Generate a LaTeX string representing the equation mathematically.
@@ -243,7 +261,7 @@ def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_ch
         str
             The LaTeX string representing the equation.
         """
-        lhs = self.lhs_term.terms[0].latex
+        lhs = self.lhs_terms[0].latex
         if drop_first_lhs_char:
             lhs = lhs[2:]
         if enclose_lhs:
@@ -295,3 +313,71 @@ def __repr__(self):
 
     def __str__(self):
         return self.__repr__()
+
+
+class ListOfAdditiveArithmeticTerms(list):
+    """Class holding list of additive arithmetic terms in equations."""
+
+    def compute_inner_products(self, basis, numerical=False, timeout=None, num_threads=None, permute=False):
+        """Compute the inner products tensor of the all the terms of the list.
+
+        Parameters
+        ----------
+        basis: SymbolicBasis
+            Basis with which to compute the inner products.
+        numerical: bool, optional
+            Whether the resulting computed inner products must be numerical or symbolic.
+            Default to `False`, i.e. symbolic output.
+        num_threads: int or None, optional
+            Number of threads to use to compute the inner products. If `None` use all the cpus available.
+            Default to `None`.
+        timeout: int or float or bool or None, optional
+            The timeout for the numerical computation of each inner product.
+            After the timeout, compute the inner product with a quadrature instead of symbolic integration.
+            Does not apply to symbolic output computations.
+            If `None` or `False`, no timeout occurs.
+            Default to `None`.
+        permute: bool, optional
+            If `True`, applies all the possible permutations to the tensor indices
+            from 1 to the rank of the tensor.
+            Default to `False`, i.e. no permutation is applied.
+        """
+
+        for term in self:
+            term.compute_inner_products(basis, numerical, timeout, num_threads, permute)
+
+    @property
+    def maximum_rank(self):
+        """int: Maximum rank of the right-hand side terms tensors."""
+        max_rank = 0
+        for term in self:
+            max_rank = max(max_rank, term.rank)
+        return max_rank
+
+    @property
+    def same_rank(self):
+        """bool: Check if all terms have the same rank."""
+        rank = self[0].rank
+        for term in self[1:]:
+            if term.rank != rank:
+                return False
+        return True
+
+    @property
+    def symbolic_expression(self):
+        """~sympy.core.expr.Expr: Symbolic expression of the collection of additive arithmetic terms."""
+
+        rterm = S.Zero
+        for term in self:
+            rterm += term.symbolic_expression
+        return rterm
+
+    @property
+    def numerical_expression(self):
+        """~sympy.core.expr.Expr: Expression of the collection of additive arithmetic terms with
+        parameters replaced by their configured values."""
+
+        rterm = S.Zero
+        for term in self:
+            rterm += term.numerical_expression
+        return rterm
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 02a07ea..52bff51 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -171,8 +171,8 @@ def compute_inner_products(self, numerical=True, timeout=None, num_threads=None)
             Default to `None`.
         """
         for field, eq in zip(self.fields, self.equations):
-            eq.lhs_term.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
-            for term in eq.terms:
+            eq.lhs_terms.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
+            for term in eq.rhs_terms:
                 term.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
 
     def compute_tensor(self, numerical=True, compute_inner_products=False, compute_inner_products_kwargs=None,
@@ -235,12 +235,14 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             for field, eq in zip(self.fields, self.equations):
                 ndim = field.state.__len__()
                 try:
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_term.inner_products.todense())
+                    if eq.other_fields_in_lhs:
+                        raise NotImplementedError('Other fields in LHS of equations is not yer implemented.')
+                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_terms.inner_products.todense())
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
-                for equation_term in eq.terms:
+                for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
-                    for term in equation_term.terms:
+                    for term in equation_term.rhs_terms:
                         term_field = term.field
                         if term_field.dynamical:
                             if self._cake is not None:
@@ -302,12 +304,12 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                             b_subs.append(sbsb)
                 ndim = field.state.__len__()
                 try:
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_term.inner_products.inverse().simplify()
+                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_terms.inner_products.inverse().simplify()
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
-                for equation_term in eq.terms:
+                for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
-                    for term in equation_term.terms:
+                    for term in equation_term.rhs_terms:
                         term_field = term.field
                         if term_field.dynamical:
                             if self._cake is not None:

From 23a766466e3cb35e6cd301d3760beb1926a93ef7 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sat, 28 Feb 2026 19:34:59 +0100
Subject: [PATCH 06/26] Adding a property to get addition of inner products for
 the LHS

---
 layercake/arithmetic/equation.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 128563f..fcf74d7 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -204,10 +204,19 @@ def numerical_lhs(self):
 
     @property
     def lhs_inner_products(self):
-        """list(~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.tensor.array.ImmutableSparseNDimArray) or list(sparse.COO(float)): Left-hand
-        side inner products of the equation, if available."""
+        """list(~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.tensor.array.ImmutableSparseNDimArray) or list(sparse.COO(float)): Inner products of each term of
+        the left-hand side of the equation, if available."""
         return [term.inner_products for term in self.lhs_terms]
 
+    @property
+    def lhs_inner_products_addition(self):
+        """~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.tensor.array.ImmutableSparseNDimArray or sparse.COO(float): Added left-hand
+        side inner products of the equation, if available. Might raise an error if not all terms are compatible."""
+        result = self.lhs_terms[0].inner_products.copy()
+        for term in self.lhs_terms[1:]:
+            result = result + term.inner_products
+        return result
+
     @property
     def maximum_rank(self):
         """int: Maximum rank of the right-hand side terms tensors."""

From 827b0adc2620225247ae15590b1be24a6706094f Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sat, 28 Feb 2026 19:41:53 +0100
Subject: [PATCH 07/26] Adapted layers to new equations org

---
 layercake/arithmetic/equation.py | 37 +++++++++++++++++++++-----------
 layercake/bakery/layers.py       | 18 +++++++++++-----
 2 files changed, 38 insertions(+), 17 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index fcf74d7..620aa64 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -222,8 +222,8 @@ def maximum_rank(self):
         """int: Maximum rank of the right-hand side terms tensors."""
         return max(self.terms.maximum_rank, self.lhs_terms.maximum_rank)
 
-    def compute_lhs_inner_products(self, basis, numerical=False, timeout=None, num_threads=None, permute=False):
-        """Compute the inner products tensor of the left-hand side terms.
+    def compute_inner_products(self, basis, numerical=False, timeout=None, num_threads=None, permute=False):
+        """Compute the inner products tensor of the left-hand and right-hand side terms.
 
         Parameters
         ----------
@@ -236,17 +236,24 @@ def compute_lhs_inner_products(self, basis, numerical=False, timeout=None, num_t
             Number of threads to use to compute the inner products. If `None` use all the cpus available.
             Default to `None`.
         timeout: int or float or bool or None, optional
-            The timeout for the numerical computation of each inner product.
-            After the timeout, compute the inner product with a quadrature instead of symbolic integration.
-            Does not apply to symbolic output computations.
-            If `None` or `False`, no timeout occurs.
-            Default to `None`.
+            Control the switch from symbolic to numerical integration. By default, `parallel_integration` workers will try to integrate
+            |Sympy| expressions symbolically, but a fallback to numerical integration can be enforced.
+            The options are:
+
+            * `None`: This is the "full-symbolic" mode. No timeout will be applied, and the switch to numerical integration will never happen.
+              Can result in very long and improbable computation time.
+            * `True`: This is the "full-numerical" mode. Symbolic computations do not occur, and the workers try directly to integrate
+              numerically.
+            * `False`: Same as `None`.
+            * An integer: defines a timeout after which, if a symbolic integration have not completed, the worker switch to the
+              numerical integration.
         permute: bool, optional
             If `True`, applies all the possible permutations to the tensor indices
             from 1 to the rank of the tensor.
             Default to `False`, i.e. no permutation is applied.
         """
         self.lhs_terms.compute_inner_products(basis, numerical, timeout, num_threads, permute)
+        self.rhs_terms.compute_inner_products(basis, numerical, timeout, num_threads, permute)
 
     def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_char=False):
         """Generate a LaTeX string representing the equation mathematically.
@@ -341,11 +348,17 @@ def compute_inner_products(self, basis, numerical=False, timeout=None, num_threa
             Number of threads to use to compute the inner products. If `None` use all the cpus available.
             Default to `None`.
         timeout: int or float or bool or None, optional
-            The timeout for the numerical computation of each inner product.
-            After the timeout, compute the inner product with a quadrature instead of symbolic integration.
-            Does not apply to symbolic output computations.
-            If `None` or `False`, no timeout occurs.
-            Default to `None`.
+            Control the switch from symbolic to numerical integration. By default, `parallel_integration` workers will try to integrate
+            |Sympy| expressions symbolically, but a fallback to numerical integration can be enforced.
+            The options are:
+
+            * `None`: This is the "full-symbolic" mode. No timeout will be applied, and the switch to numerical integration will never happen.
+              Can result in very long and improbable computation time.
+            * `True`: This is the "full-numerical" mode. Symbolic computations do not occur, and the workers try directly to integrate
+              numerically.
+            * `False`: Same as `None`.
+            * An integer: defines a timeout after which, if a symbolic integration have not completed, the worker switch to the
+              numerical integration.
         permute: bool, optional
             If `True`, applies all the possible permutations to the tensor indices
             from 1 to the rank of the tensor.
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 52bff51..dc794d0 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -156,7 +156,7 @@ def maximum_rank(self):
 
     def compute_inner_products(self, numerical=True, timeout=None, num_threads=None):
         """Compute the inner products tensors, either symbolic or numerical ones, of all the terms
-        of the layer equations, including the left-hand side term.
+        of the layer equations, including the left-hand side terms.
         Computations are parallelized on multiple CPUs.
 
         Parameters
@@ -165,15 +165,23 @@ def compute_inner_products(self, numerical=True, timeout=None, num_threads=None)
             Whether to compute numerical or symbolic inner products.
             Default to `True` (numerical inner products as output).
         timeout: int or bool or None, optional
-            TODO
+            Control the switch from symbolic to numerical integration. By default, `parallel_integration` workers will try to integrate
+            |Sympy| expressions symbolically, but a fallback to numerical integration can be enforced.
+            The options are:
+
+            * `None`: This is the "full-symbolic" mode. No timeout will be applied, and the switch to numerical integration will never happen.
+              Can result in very long and improbable computation time.
+            * `True`: This is the "full-numerical" mode. Symbolic computations do not occur, and the workers try directly to integrate
+              numerically.
+            * `False`: Same as `None`.
+            * An integer: defines a timeout after which, if a symbolic integration have not completed, the worker switch to the
+              numerical integration.
         num_threads: None or int, optional
             Number of CPUs to use in parallel for the computations. If `None`, use all the CPUs available.
             Default to `None`.
         """
         for field, eq in zip(self.fields, self.equations):
-            eq.lhs_terms.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
-            for term in eq.rhs_terms:
-                term.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
+            eq.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
 
     def compute_tensor(self, numerical=True, compute_inner_products=False, compute_inner_products_kwargs=None,
                        substitutions=None, basis_subs=False, parameters_subs=None):

From 13f6f2488cb1bb4a6f42df5779e5f1ddf36ed78e Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sat, 28 Feb 2026 19:48:03 +0100
Subject: [PATCH 08/26] Cont'd

---
 examples/atmospheric/baroclinic_one_layer.py         | 4 ++--
 examples/atmospheric/barotropic_one_layer.py         | 2 +-
 examples/atmospheric/stratospheric_planetary_flow.py | 2 +-
 examples/ocean-atmosphere/maooam.py                  | 8 ++++----
 layercake/bakery/layers.py                           | 6 ++++--
 test/test_RP1982_qgs_numerical.py                    | 4 ++--
 test/test_RP1982_qgs_symbolic.py                     | 4 ++--
 test/test_maooam_qgs_numerical.py                    | 8 ++++----
 test/test_maooam_qgs_symbolic.py                     | 8 ++++----
 9 files changed, 24 insertions(+), 22 deletions(-)

diff --git a/examples/atmospheric/baroclinic_one_layer.py b/examples/atmospheric/baroclinic_one_layer.py
index 9c327d8..b85ce41 100644
--- a/examples/atmospheric/baroclinic_one_layer.py
+++ b/examples/atmospheric/baroclinic_one_layer.py
@@ -122,7 +122,7 @@
 
 # defining the LHS as the time derivative of the vorticity
 vorticity = OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system)
-barotropic_equation = Equation(psi, lhs_term=vorticity)
+barotropic_equation = Equation(psi, lhs_terms=vorticity)
 
 # Defining the advection term
 advection_term1 = vorticity_advection(psi, psi, atmospheric_basis.coordinate_system, sign=-1)
@@ -165,7 +165,7 @@
 
 lin_lhs = LinearTerm(theta, prefactor=a, sign=-1)
 lhs = AdditionOfTerms(lin_lhs, vorticity)
-baroclinic_equation = Equation(theta, lhs_term=lhs)
+baroclinic_equation = Equation(theta, lhs_terms=lhs)
 
 # Defining the advection terms
 advection_term1 = vorticity_advection(psi, theta, atmospheric_basis.coordinate_system, sign=-1)
diff --git a/examples/atmospheric/barotropic_one_layer.py b/examples/atmospheric/barotropic_one_layer.py
index 70c8cfc..8d2c52e 100644
--- a/examples/atmospheric/barotropic_one_layer.py
+++ b/examples/atmospheric/barotropic_one_layer.py
@@ -75,7 +75,7 @@
 
 # defining the LHS as the time derivative of the vorticity
 vorticity = OperatorTerm(psi, Laplacian, basis.coordinate_system)
-barotropic_equation = Equation(psi, lhs_term=vorticity)
+barotropic_equation = Equation(psi, lhs_terms=vorticity)
 
 # defining the advection term
 advection_term = vorticity_advection(psi, psi, basis.coordinate_system, sign=-1)
diff --git a/examples/atmospheric/stratospheric_planetary_flow.py b/examples/atmospheric/stratospheric_planetary_flow.py
index e702401..da7bec6 100644
--- a/examples/atmospheric/stratospheric_planetary_flow.py
+++ b/examples/atmospheric/stratospheric_planetary_flow.py
@@ -61,7 +61,7 @@
 
 # defining the LHS as the time derivative of the vorticity
 vorticity = OperatorTerm(psi, Laplacian, basis.coordinate_system)
-planetary_equation = Equation(psi, lhs_term=vorticity)
+planetary_equation = Equation(psi, lhs_terms=vorticity)
 
 # defining the advection term
 advection_term = vorticity_advection(psi, psi, basis.coordinate_system, sign=-1)
diff --git a/examples/ocean-atmosphere/maooam.py b/examples/ocean-atmosphere/maooam.py
index 4145691..5352973 100644
--- a/examples/ocean-atmosphere/maooam.py
+++ b/examples/ocean-atmosphere/maooam.py
@@ -155,7 +155,7 @@
 
 # defining LHS as the time derivative of the vorticity
 vorticity = OperatorTerm(psi_a, Laplacian, atmospheric_basis.coordinate_system)
-barotropic_equation = Equation(psi_a, lhs_term=vorticity)
+barotropic_equation = Equation(psi_a, lhs_terms=vorticity)
 
 # Defining the advection term
 advection_term1 = vorticity_advection(psi_a, psi_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -190,7 +190,7 @@
 
 lin_lhs = LinearTerm(theta_a, prefactor=a, sign=-1)
 lhs = AdditionOfTerms(lin_lhs, vorticity)
-baroclinic_equation = Equation(theta_a, lhs_term=lhs)
+baroclinic_equation = Equation(theta_a, lhs_terms=lhs)
 
 # Defining the advection terms
 advection_term1 = vorticity_advection(psi_a, theta_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -250,7 +250,7 @@
 vorticity = OperatorTerm(psi_o, Laplacian, oceanic_basis.coordinate_system)
 lin_lhs = LinearTerm(psi_o, prefactor=G)
 lhs = AdditionOfTerms(lin_lhs, vorticity)
-oceanic_equation = Equation(psi_o, lhs_term=lhs)
+oceanic_equation = Equation(psi_o, lhs_terms=lhs)
 
 # Defining the advection term
 advection_term = vorticity_advection(psi_o, psi_o, oceanic_basis.coordinate_system, sign=-1)
@@ -283,7 +283,7 @@
 
 # Defining the LHS
 lhs = LinearTerm(deltaT_o)
-ocean_temperature_equation = Equation(deltaT_o, lhs_term=lhs)
+ocean_temperature_equation = Equation(deltaT_o, lhs_terms=lhs)
 
 # Defining the advection term
 advection_term = Jacobian(psi_o, deltaT_o, oceanic_basis.coordinate_system, sign=-1)
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index dc794d0..8b79512 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -245,7 +245,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 try:
                     if eq.other_fields_in_lhs:
                         raise NotImplementedError('Other fields in LHS of equations is not yer implemented.')
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_terms.inner_products.todense())
+                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
@@ -312,7 +312,9 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                             b_subs.append(sbsb)
                 ndim = field.state.__len__()
                 try:
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_terms.inner_products.inverse().simplify()
+                    if eq.other_fields_in_lhs:
+                        raise NotImplementedError('Other fields in LHS of equations is not yer implemented.')
+                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inverse().simplify()
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
diff --git a/test/test_RP1982_qgs_numerical.py b/test/test_RP1982_qgs_numerical.py
index 00f679e..5b5d5f1 100644
--- a/test/test_RP1982_qgs_numerical.py
+++ b/test/test_RP1982_qgs_numerical.py
@@ -184,7 +184,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         vorticity = OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system)
-        barotropic_equation = Equation(psi, lhs_term=vorticity)
+        barotropic_equation = Equation(psi, lhs_terms=vorticity)
 
         # Defining the advection term
         advection_term1 = vorticity_advection(psi, psi, atmospheric_basis.coordinate_system, sign=-1)
@@ -228,7 +228,7 @@ def layercake_outputs(self, output_func=None):
 
         lin_lhs = LinearTerm(theta, prefactor=a, sign=-1)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        baroclinic_equation = Equation(theta, lhs_term=lhs)
+        baroclinic_equation = Equation(theta, lhs_terms=lhs)
 
         # Defining the advection terms
         advection_term1 = vorticity_advection(psi, theta, atmospheric_basis.coordinate_system, sign=-1)
diff --git a/test/test_RP1982_qgs_symbolic.py b/test/test_RP1982_qgs_symbolic.py
index 46fd544..5e66354 100644
--- a/test/test_RP1982_qgs_symbolic.py
+++ b/test/test_RP1982_qgs_symbolic.py
@@ -185,7 +185,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         vorticity = OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system)
-        barotropic_equation = Equation(psi, lhs_term=vorticity)
+        barotropic_equation = Equation(psi, lhs_terms=vorticity)
 
         # Defining the advection term
         advection_term1 = vorticity_advection(psi, psi, atmospheric_basis.coordinate_system, sign=-1)
@@ -229,7 +229,7 @@ def layercake_outputs(self, output_func=None):
 
         lin_lhs = LinearTerm(theta, prefactor=a, sign=-1)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        baroclinic_equation = Equation(theta, lhs_term=lhs)
+        baroclinic_equation = Equation(theta, lhs_terms=lhs)
 
         # Defining the advection terms
         advection_term1 = vorticity_advection(psi, theta, atmospheric_basis.coordinate_system, sign=-1)
diff --git a/test/test_maooam_qgs_numerical.py b/test/test_maooam_qgs_numerical.py
index 821e9a4..efe4ade 100644
--- a/test/test_maooam_qgs_numerical.py
+++ b/test/test_maooam_qgs_numerical.py
@@ -227,7 +227,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         vorticity = OperatorTerm(psi_a, Laplacian, atmospheric_basis.coordinate_system)
-        barotropic_equation = Equation(psi_a, lhs_term=vorticity)
+        barotropic_equation = Equation(psi_a, lhs_terms=vorticity)
 
         # Defining the advection term
         advection_term1 = vorticity_advection(psi_a, psi_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -262,7 +262,7 @@ def layercake_outputs(self, output_func=None):
 
         lin_lhs = LinearTerm(theta_a, prefactor=a, sign=-1)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        baroclinic_equation = Equation(theta_a, lhs_term=lhs)
+        baroclinic_equation = Equation(theta_a, lhs_terms=lhs)
 
         # Defining the advection terms
         advection_term1 = vorticity_advection(psi_a, theta_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -322,7 +322,7 @@ def layercake_outputs(self, output_func=None):
         vorticity = OperatorTerm(psi_o, Laplacian, oceanic_basis.coordinate_system)
         lin_lhs = LinearTerm(psi_o, prefactor=G)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        oceanic_equation = Equation(psi_o, lhs_term=lhs)
+        oceanic_equation = Equation(psi_o, lhs_terms=lhs)
 
         # Defining the advection term
         advection_term = vorticity_advection(psi_o, psi_o, oceanic_basis.coordinate_system, sign=-1)
@@ -355,7 +355,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         lhs = LinearTerm(deltaT_o)
-        ocean_temperature_equation = Equation(deltaT_o, lhs_term=lhs)
+        ocean_temperature_equation = Equation(deltaT_o, lhs_terms=lhs)
 
         # Defining the advection term
         advection_term = Jacobian(psi_o, deltaT_o, oceanic_basis.coordinate_system, sign=-1)
diff --git a/test/test_maooam_qgs_symbolic.py b/test/test_maooam_qgs_symbolic.py
index 515f19f..817e8b8 100644
--- a/test/test_maooam_qgs_symbolic.py
+++ b/test/test_maooam_qgs_symbolic.py
@@ -229,7 +229,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         vorticity = OperatorTerm(psi_a, Laplacian, atmospheric_basis.coordinate_system)
-        barotropic_equation = Equation(psi_a, lhs_term=vorticity)
+        barotropic_equation = Equation(psi_a, lhs_terms=vorticity)
 
         # Defining the advection term
         advection_term1 = vorticity_advection(psi_a, psi_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -264,7 +264,7 @@ def layercake_outputs(self, output_func=None):
 
         lin_lhs = LinearTerm(theta_a, prefactor=a, sign=-1)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        baroclinic_equation = Equation(theta_a, lhs_term=lhs)
+        baroclinic_equation = Equation(theta_a, lhs_terms=lhs)
 
         # Defining the advection terms
         advection_term1 = vorticity_advection(psi_a, theta_a, atmospheric_basis.coordinate_system, sign=-1)
@@ -324,7 +324,7 @@ def layercake_outputs(self, output_func=None):
         vorticity = OperatorTerm(psi_o, Laplacian, oceanic_basis.coordinate_system)
         lin_lhs = LinearTerm(psi_o, prefactor=G)
         lhs = AdditionOfTerms(lin_lhs, vorticity)
-        oceanic_equation = Equation(psi_o, lhs_term=lhs)
+        oceanic_equation = Equation(psi_o, lhs_terms=lhs)
 
         # Defining the advection term
         advection_term = vorticity_advection(psi_o, psi_o, oceanic_basis.coordinate_system, sign=-1)
@@ -357,7 +357,7 @@ def layercake_outputs(self, output_func=None):
         # Defining the equation and LHS
         # Laplacian
         lhs = LinearTerm(deltaT_o)
-        ocean_temperature_equation = Equation(deltaT_o, lhs_term=lhs)
+        ocean_temperature_equation = Equation(deltaT_o, lhs_terms=lhs)
 
         # Defining the advection term
         advection_term = Jacobian(psi_o, deltaT_o, oceanic_basis.coordinate_system, sign=-1)

From b5b914a7adecca2db81484013774a4bd4a90dd12 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 09:06:48 +0100
Subject: [PATCH 09/26] Correcting a recursion bug

---
 layercake/arithmetic/equation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 620aa64..69feb06 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -113,7 +113,7 @@ def add_lhs_terms(self, terms):
     @property
     def terms(self):
         """Alias for the list of RHS arithmetic terms."""
-        return self.terms
+        return self.rhs_terms
 
     @property
     def other_fields(self):

From b34c1723d521171e1b3fa81b57477f8cb0162d66 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 09:07:57 +0100
Subject: [PATCH 10/26] Comment

---
 layercake/arithmetic/equation.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 69feb06..4266b5e 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -128,6 +128,7 @@ def other_fields(self):
 
     @property
     def other_fields_in_lhs(self):
+        """list(~field.Field): List of additional fields present in the LHS of the equation."""
         other_fields = list()
         for term in self.lhs_terms:
             if term.field is not self.field and term.field.dynamical and term.field not in other_fields:

From f145e970f7e77ad5c11064c645985ead7b7615f4 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 09:12:01 +0100
Subject: [PATCH 11/26] Correcting previous refactoring

---
 layercake/bakery/layers.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 8b79512..b9c8104 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -250,7 +250,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
-                    for term in equation_term.rhs_terms:
+                    for term in equation_term.terms:
                         term_field = term.field
                         if term_field.dynamical:
                             if self._cake is not None:
@@ -319,7 +319,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
-                    for term in equation_term.rhs_terms:
+                    for term in equation_term.terms:
                         term_field = term.field
                         if term_field.dynamical:
                             if self._cake is not None:

From 07df4378f3871eb6e33fb31840aa60297b6f3076 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 09:17:37 +0100
Subject: [PATCH 12/26] Cont'd

---
 layercake/arithmetic/equation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 4266b5e..cd978e5 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -120,7 +120,7 @@ def other_fields(self):
         """list(~field.Field): List of additional fields present in the equation."""
         other_fields = list()
         for equation_term in self.terms:
-            for term in equation_term.rhs_terms:
+            for term in equation_term.terms:
                 if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
                     other_fields.append(term.field)
         other_fields = other_fields + self.other_fields_in_lhs
@@ -140,7 +140,7 @@ def parameter_fields(self):
         """list(~field.ParameterField): List of non-dynamical parameter fields present in the equation."""
         parameter_fields = list()
         for equation_term in self.terms:
-            for term in equation_term.rhs_terms:
+            for term in equation_term.terms:
                 if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
                     parameter_fields.append(term.field)
         for term in self.lhs_terms:

From 357ef7e4910df016624308bd708df20d02ee2912 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 11:44:51 +0100
Subject: [PATCH 13/26] Added possibility to have other fields from the same
 layer in the LHS of layer's equations.

---
 layercake/bakery/layers.py | 68 +++++++++++++++++++++++++++++++++-----
 1 file changed, 59 insertions(+), 9 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index b9c8104..cf180f4 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -37,7 +37,7 @@ class Layer(object):
     Parameters
     ----------
     name: str, optional
-        Optional name for the layer.
+        Name for the layer.
 
     Attributes
     ----------
@@ -57,6 +57,9 @@ def __init__(self, name=''):
         self.name = name
         self._cake = None
         self._cake_order = 0
+        self._lhs_inversion = True
+        self._lhs_inverted = False
+        self._lhs_mat = None
 
     @property
     def _cake_first_index(self):
@@ -105,6 +108,18 @@ def other_fields(self):
                     other_fields.append(field)
         return other_fields
 
+    @property
+    def other_fields_in_lhs(self):
+        """list(~field.Field): Returns the list of dynamical fields of other layers present in the LHS of the layer's equations,
+        i.e. the fields whose time evolution is provided by the partial differential equations of other layers."""
+        layer_fields = self.fields
+        other_fields = list()
+        for eq in self.equations:
+            for field in eq.other_fields_in_lhs:
+                if field not in layer_fields:
+                    other_fields.append(field)
+        return other_fields
+
     @property
     def parameters(self):
         """list(~parameter.Parameter): Returns the list of parameters of the layer, i.e. the explicit parameters
@@ -238,14 +253,28 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
         if numerical:
 
             self.tensor = sp.zeros(shape, dtype=np.float64, format='dok')
-            lhs_mat = np.zeros((self.ndim+1, self.ndim+1))
+            self._lhs_mat = sp.zeros((self.ndim+1, self.ndim+1), dtype=np.float64, format='dok')
+            lhs_mat_inverted = np.zeros((self.ndim+1, self.ndim+1))
             lhs_order = 1
             for field, eq in zip(self.fields, self.equations):
                 ndim = field.state.__len__()
                 try:
                     if eq.other_fields_in_lhs:
-                        raise NotImplementedError('Other fields in LHS of equations is not yer implemented.')
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
+                        if self.other_fields_in_lhs:
+                            self._lhs_inversion = False
+                            raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
+                        else:
+                            for lhs_term in eq.lhs_terms:
+                                ofield = lhs_term.field
+                                ofield_order = 1
+                                for tfield in self.fields:
+                                    if ofield is tfield:
+                                        break
+                                    ofield_order += ndim
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products
+                    else:
+                        lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
+                        self._lhs_inverted = True
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
@@ -288,7 +317,10 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     if np.any(increment != 0):
                         self.tensor[args] = self.tensor[args] + increment
                 lhs_order += ndim
-            self.tensor = sp.COO(np.tensordot(lhs_mat, self.tensor.to_coo(), 1))
+            if self._lhs_inversion and not self._lhs_inverted:
+                lhs_mat_inverted = np.linalg.inv(self._lhs_mat.todense())
+                self._lhs_inverted = True
+            self.tensor = sp.COO(np.tensordot(lhs_mat_inverted, self.tensor.to_coo(), 1))
 
         else:
             b_subs = list()
@@ -299,7 +331,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             else:
                 p_subs = list()
             self.tensor = MutableSparseNDimArray(iterable={}, shape=shape)
-            lhs_mat = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
+            self._lhs_mat = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
+            lhs_mat_inverted = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
             lhs_order = 1
             for field, eq in zip(self.fields, self.equations):
                 bsb = field.basis.substitutions
@@ -313,8 +346,21 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 ndim = field.state.__len__()
                 try:
                     if eq.other_fields_in_lhs:
-                        raise NotImplementedError('Other fields in LHS of equations is not yer implemented.')
-                    lhs_mat[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inverse().simplify()
+                        if self.other_fields_in_lhs:
+                            self._lhs_inversion = False
+                            raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
+                        else:
+                            for lhs_term in eq.lhs_terms:
+                                ofield = lhs_term.field
+                                ofield_order = 1
+                                for tfield in self.fields:
+                                    if ofield is tfield:
+                                        break
+                                    ofield_order += ndim
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products
+                    else:
+                        lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inv().simplify()
+                        self._lhs_inverted = True
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
                 for equation_term in eq.rhs_terms:
@@ -380,7 +426,11 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                         increment = ImmutableSparseNDimArray(increment)
                     self.tensor[args] = self.tensor[args] + increment
                 lhs_order += ndim
-            self.tensor = (ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat, self.tensor, 1))
+
+            if self._lhs_inversion and not self._lhs_inverted:
+                lhs_mat_inverted = self._lhs_mat.inv().simplify()
+                self._lhs_inverted = True
+            self.tensor = (ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat_inverted, self.tensor, 1))
                            .subs(b_subs).subs(p_subs).subs(substitutions))
 
     def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_char=False):

From b8df6546ca5cbdc3baee57fc1c4d3b4b6a944136 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 11:53:02 +0100
Subject: [PATCH 14/26] Corrected a bug in new other fields in lhs property

---
 layercake/arithmetic/equation.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index cd978e5..30a74a1 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -130,9 +130,10 @@ def other_fields(self):
     def other_fields_in_lhs(self):
         """list(~field.Field): List of additional fields present in the LHS of the equation."""
         other_fields = list()
-        for term in self.lhs_terms:
-            if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
-                other_fields.append(term.field)
+        for equation_term in self.lhs_terms:
+            for term in equation_term.terms:
+                if term.field is not self.field and term.field.dynamical and term.field not in other_fields:
+                    other_fields.append(term.field)
         return other_fields
 
     @property

From 20962d8fd77a018ff7d86615b960978d3742d156 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 1 Mar 2026 12:04:22 +0100
Subject: [PATCH 15/26] Corrected same bug in parameter_fields property

---
 layercake/arithmetic/equation.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/layercake/arithmetic/equation.py b/layercake/arithmetic/equation.py
index 30a74a1..4c17e0e 100644
--- a/layercake/arithmetic/equation.py
+++ b/layercake/arithmetic/equation.py
@@ -144,9 +144,10 @@ def parameter_fields(self):
             for term in equation_term.terms:
                 if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
                     parameter_fields.append(term.field)
-        for term in self.lhs_terms:
-            if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
-                parameter_fields.append(term.field)
+        for equation_term in self.lhs_terms:
+            for term in equation_term.terms:
+                if term.field is not self.field and not term.field.dynamical and term.field not in parameter_fields:
+                    parameter_fields.append(term.field)
         return parameter_fields
 
     @property

From a620b244b87b3357f109c7e27b3ab64152fb7066 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Sun, 8 Mar 2026 21:51:12 +0100
Subject: [PATCH 16/26] Corrected inversion of layers LHS issues

---
 layercake/bakery/layers.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index cf180f4..b9dc8dd 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -271,7 +271,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                     if ofield is tfield:
                                         break
                                     ofield_order += ndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products.todense()
                     else:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
                         self._lhs_inverted = True
@@ -318,7 +318,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                         self.tensor[args] = self.tensor[args] + increment
                 lhs_order += ndim
             if self._lhs_inversion and not self._lhs_inverted:
-                lhs_mat_inverted = np.linalg.inv(self._lhs_mat.todense())
+                lhs_mat_inverted[1:, 1:] = np.linalg.inv(self._lhs_mat.todense()[1:, 1:])
                 self._lhs_inverted = True
             self.tensor = sp.COO(np.tensordot(lhs_mat_inverted, self.tensor.to_coo(), 1))
 

From 7fdd2f1ee2c1c7999c00520f19542d11cbbf6846 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Mon, 9 Mar 2026 13:58:34 +0100
Subject: [PATCH 17/26] Cont'd

---
 layercake/bakery/layers.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index b9dc8dd..9cbd1df 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -271,7 +271,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                     if ofield is tfield:
                                         break
                                     ofield_order += ndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products.todense()
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = \
+                                    self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] + lhs_term.inner_products.todense()
                     else:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
                         self._lhs_inverted = True
@@ -357,7 +358,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                     if ofield is tfield:
                                         break
                                     ofield_order += ndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = lhs_term.inner_products
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] += lhs_term.inner_products
                     else:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inv().simplify()
                         self._lhs_inverted = True

From 51979ddfbae153ec7a25ac1ab16d30b6492285f8 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 09:58:06 +0100
Subject: [PATCH 18/26] Error management

---
 layercake/bakery/layers.py | 83 +++++++++++++++++++++-----------------
 1 file changed, 46 insertions(+), 37 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 9cbd1df..c249a0d 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -18,6 +18,7 @@
 import matplotlib.pyplot as plt
 import sparse as sp
 from copy import deepcopy
+import warnings
 
 from sympy import MutableSparseNDimArray, MutableSparseMatrix, ImmutableMatrix, ImmutableSparseNDimArray
 from sympy import zeros as sympy_zeros
@@ -258,26 +259,27 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             lhs_order = 1
             for field, eq in zip(self.fields, self.equations):
                 ndim = field.state.__len__()
-                try:
-                    if eq.other_fields_in_lhs:
-                        if self.other_fields_in_lhs:
-                            self._lhs_inversion = False
-                            raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
-                        else:
-                            for lhs_term in eq.lhs_terms:
-                                ofield = lhs_term.field
-                                ofield_order = 1
-                                for tfield in self.fields:
-                                    if ofield is tfield:
-                                        break
-                                    ofield_order += ndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = \
-                                    self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] + lhs_term.inner_products.todense()
+                if eq.other_fields_in_lhs:
+                    if self.other_fields_in_lhs:
+                        self._lhs_inversion = False
+                        raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
                     else:
+                        for lhs_term in eq.lhs_terms:
+                            ofield = lhs_term.field
+                            ofield_order = 1
+                            for tfield in self.fields:
+                                if ofield is tfield:
+                                    break
+                                ofield_order += ndim
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = \
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] + lhs_term.inner_products.todense()
+                else:
+                    try:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
                         self._lhs_inverted = True
-                except LinAlgError:
-                    raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
+                    except LinAlgError:
+                        raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
+
                 for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
                     for term in equation_term.terms:
@@ -319,8 +321,11 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                         self.tensor[args] = self.tensor[args] + increment
                 lhs_order += ndim
             if self._lhs_inversion and not self._lhs_inverted:
-                lhs_mat_inverted[1:, 1:] = np.linalg.inv(self._lhs_mat.todense()[1:, 1:])
-                self._lhs_inverted = True
+                try:
+                    lhs_mat_inverted[1:, 1:] = np.linalg.inv(self._lhs_mat.todense()[1:, 1:])
+                    self._lhs_inverted = True
+                except LinAlgError:
+                    raise LinAlgError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')
             self.tensor = sp.COO(np.tensordot(lhs_mat_inverted, self.tensor.to_coo(), 1))
 
         else:
@@ -345,25 +350,26 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                         else:
                             b_subs.append(sbsb)
                 ndim = field.state.__len__()
-                try:
-                    if eq.other_fields_in_lhs:
-                        if self.other_fields_in_lhs:
-                            self._lhs_inversion = False
-                            raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
-                        else:
-                            for lhs_term in eq.lhs_terms:
-                                ofield = lhs_term.field
-                                ofield_order = 1
-                                for tfield in self.fields:
-                                    if ofield is tfield:
-                                        break
-                                    ofield_order += ndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] += lhs_term.inner_products
+                if eq.other_fields_in_lhs:
+                    if self.other_fields_in_lhs:
+                        self._lhs_inversion = False
+                        raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
                     else:
+                        for lhs_term in eq.lhs_terms:
+                            ofield = lhs_term.field
+                            ofield_order = 1
+                            for tfield in self.fields:
+                                if ofield is tfield:
+                                    break
+                                ofield_order += ndim
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] += lhs_term.inner_products
+                else:
+                    try:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inv().simplify()
                         self._lhs_inverted = True
-                except NonInvertibleMatrixError:
-                    raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
+                    except NonInvertibleMatrixError:
+                        raise NonInvertibleMatrixError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
+
                 for equation_term in eq.rhs_terms:
                     slices = [slice(lhs_order, lhs_order + ndim)]
                     for term in equation_term.terms:
@@ -429,8 +435,11 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 lhs_order += ndim
 
             if self._lhs_inversion and not self._lhs_inverted:
-                lhs_mat_inverted = self._lhs_mat.inv().simplify()
-                self._lhs_inverted = True
+                try:
+                    lhs_mat_inverted = self._lhs_mat.inv().simplify()
+                    self._lhs_inverted = True
+                except NonInvertibleMatrixError:
+                    raise NonInvertibleMatrixError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')
             self.tensor = (ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat_inverted, self.tensor, 1))
                            .subs(b_subs).subs(p_subs).subs(substitutions))
 

From 381efb0dcba9e58634ebbb5a339c17d29b5468b0 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 10:23:07 +0100
Subject: [PATCH 19/26] Added case when LHS inversion must be cake-wise

---
 layercake/bakery/layers.py | 69 ++++++++++++++++++++++++++++++++------
 1 file changed, 59 insertions(+), 10 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index c249a0d..af82c0d 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -262,17 +262,38 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 if eq.other_fields_in_lhs:
                     if self.other_fields_in_lhs:
                         self._lhs_inversion = False
-                        raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
+                        if self._cake is None:
+                            raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
+                        else:
+                            if np.all(self._lhs_mat.todense() == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                                self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
+                            if np.any(self._lhs_mat.todense() != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                                raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+
+                        for lhs_term in eq.lhs_terms:
+                            ofield = lhs_term.field
+                            ofield_order = 1
+                            ondim = ofield.state.__len__()
+                            for tfield in self._cake.fields:
+                                if ofield is tfield:
+                                    break
+                                tndim = tfield.state.__len__()
+                                ofield_order += tndim
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] = \
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
+
                     else:
                         for lhs_term in eq.lhs_terms:
                             ofield = lhs_term.field
                             ofield_order = 1
+                            ondim = ofield.state.__len__()
                             for tfield in self.fields:
                                 if ofield is tfield:
                                     break
-                                ofield_order += ndim
-                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] = \
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] + lhs_term.inner_products.todense()
+                                tndim = tfield.state.__len__()
+                                ofield_order += tndim
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] = \
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
                 else:
                     try:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
@@ -326,7 +347,11 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     self._lhs_inverted = True
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')
-            self.tensor = sp.COO(np.tensordot(lhs_mat_inverted, self.tensor.to_coo(), 1))
+            if self._lhs_inverted:
+                self.tensor = sp.COO(np.tensordot(lhs_mat_inverted, self.tensor.to_coo(), 1))
+            else:
+                warnings.warn(f'LHS of layer {self} has not been inverted. Its tensor hold the RHS tendencies alone.')
+                self.tensor = self.tensor.to_coo()
 
         else:
             b_subs = list()
@@ -353,16 +378,35 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 if eq.other_fields_in_lhs:
                     if self.other_fields_in_lhs:
                         self._lhs_inversion = False
-                        raise NotImplementedError('Fields from other layers in LHS of equations is not yer implemented.')
+                        if self._cake is None:
+                            raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
+                        else:
+                            if np.all(self._lhs_mat == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                                self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
+                            if np.any(self._lhs_mat != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                                raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+
+                        for lhs_term in eq.lhs_terms:
+                            ofield = lhs_term.field
+                            ofield_order = 1
+                            ondim = ofield.state.__len__()
+                            for tfield in self._cake.fields:
+                                if ofield is tfield:
+                                    break
+                                tndim = tfield.state.__len__()
+                                ofield_order += tndim
+                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
                     else:
                         for lhs_term in eq.lhs_terms:
                             ofield = lhs_term.field
                             ofield_order = 1
+                            ondim = ofield.state.__len__()
                             for tfield in self.fields:
                                 if ofield is tfield:
                                     break
-                                ofield_order += ndim
-                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ndim] += lhs_term.inner_products
+                                tndim = tfield.state.__len__()
+                                ofield_order += tndim
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
                 else:
                     try:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inv().simplify()
@@ -440,8 +484,13 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     self._lhs_inverted = True
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')
-            self.tensor = (ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat_inverted, self.tensor, 1))
-                           .subs(b_subs).subs(p_subs).subs(substitutions))
+            if self._lhs_inverted:
+                self.tensor = (ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat_inverted, self.tensor, 1))
+                               .subs(b_subs).subs(p_subs).subs(substitutions))
+            else:
+                warnings.warn(f'LHS of layer {self} has not been inverted. Its tensor hold the RHS tendencies alone.')
+                self.tensor = ImmutableSparseNDimArray(self.tensor).subs(b_subs).subs(p_subs).subs(substitutions)
+                self._lhs_mat = self._lhs_mat.subs(b_subs).subs(p_subs).subs(substitutions)
 
     def to_latex(self, enclose_lhs=True, drop_first_lhs_char=True, drop_first_rhs_char=False):
         """Generate the LaTeX strings representing the layer's equations mathematically.

From c7c9d7771f50410cdf1bff1eecd6aff1b920bd1f Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 10:50:53 +0100
Subject: [PATCH 20/26] Cont'd (to be tested)

---
 layercake/bakery/cake.py   | 46 +++++++++++++++++++++++++++++++++-----
 layercake/bakery/layers.py | 10 +++++----
 2 files changed, 46 insertions(+), 10 deletions(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index 5ec9cfe..4eb6da2 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -15,18 +15,25 @@
 
 from contextlib import redirect_stdout
 import numpy as np
+from numpy.linalg import LinAlgError
 import matplotlib.pyplot as plt
 import sparse as sp
 from numba import njit
+
+from sympy import ImmutableSparseNDimArray, MutableSparseNDimArray
+from sympy import MutableSparseMatrix
+from sympy import zeros as sympy_zeros
+from sympy import simplify, N
+from sympy.tensor.array import permutedims
+from sympy.matrices.exceptions import NonInvertibleMatrixError
+
 from layercake.utils.tensor import sparse_mul, jsparse_mul
 from layercake.utils.symbolic_tensor import get_coords_and_values_from_tensor, compute_jacobian_permutations
 from layercake.formatters.fortran import FortranJacobianEquationFormatter, FortranEquationFormatter
 from layercake.formatters.python import PythonJacobianEquationFormatter, PythonEquationFormatter
 from layercake.formatters.julia import JuliaJacobianEquationFormatter, JuliaEquationFormatter
 from layercake.utils import isin
-from sympy import ImmutableSparseNDimArray, MutableSparseNDimArray
-from sympy import simplify, N
-from sympy.tensor.array import permutedims
+from layercake.utils.symbolic_tensor import symbolic_tensordot
 
 real_eps = np.finfo(np.float64).eps
 small_number = 1.e-12
@@ -45,6 +52,7 @@ class Cake(object):
     def __init__(self):
 
         self.layers = list()
+        self._lhs_inversion = True
 
     def add_layer(self, layer):
         """Add a layer object to the cake.
@@ -150,6 +158,11 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             Only applies for the symbolic tendencies.
 
         """
+        self._lhs_inversion = True
+        for layer in self.layers:
+            if not layer._lhs_inversion:
+                self._lhs_inversion = False
+                break
 
         for layer in self.layers:
             layer.compute_tensor(numerical=numerical,
@@ -157,7 +170,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                  compute_inner_products_kwargs=compute_inner_products_kwargs,
                                  substitutions=substitutions,
                                  basis_subs=basis_subs,
-                                 parameters_subs=parameters_subs
+                                 parameters_subs=parameters_subs,
+                                 lhs_inversion=self._lhs_inversion
                                  )
 
     @property
@@ -205,8 +219,12 @@ def tensor(self):
 
         if numerical:
             tensor = sp.zeros(shape, dtype=np.float64, format='dok')
+            if self._lhs_inversion:
+                lhs_mat = sp.zeros((self.ndim + 1, self.ndim + 1), dtype=np.float64, format='dok')
         else:
             tensor = MutableSparseNDimArray(iterable={}, shape=shape)
+            if self._lhs_inversion:
+                lhs_mat = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
 
         for i, layer in enumerate(self.layers):
             if (numerical and not isinstance(layer.tensor, sp.COO) or
@@ -225,13 +243,29 @@ def tensor(self):
             args = tuple(slices + zeros)
             if numerical:
                 tensor[args] = tensor[args] + layer.tensor.todense()[1:]
+                if self._lhs_inversion:
+                    lhs_mat[args[:2]] = lhs_mat[args[:1]] + layer._lhs_mat.todense()[1:]
             else:
                 tensor[args] = tensor[args] + layer.tensor[1:]
+                if self._lhs_inversion:
+                    lhs_mat[args[:2]] = lhs_mat[args[:2]] + layer._lhs_mat[1:]
 
         if numerical:
-            tensor = tensor.to_coo()
+            if self._lhs_inversion:
+                try:
+                    tensor = sp.COO(np.tensordot(np.linalg.inv(lhs_mat), tensor.to_coo(), 1))
+                except LinAlgError:
+                    raise LinAlgError(f'The left-hand side of the cake is not invertible with the provided basis.')
+            else:
+                tensor = tensor.to_coo()
         else:
-            tensor = ImmutableSparseNDimArray(tensor)
+            if self._lhs_inversion:
+                try:
+                    tensor = ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat.inv().simplify(), tensor, 1))
+                except NonInvertibleMatrixError:
+                    raise NonInvertibleMatrixError(f'The left-hand side of the cake is not invertible with the provided basis.')
+            else:
+                tensor = ImmutableSparseNDimArray(tensor)
 
         tensor = self.simplify_tensor(tensor, small_number)
 
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index af82c0d..39aa00d 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -200,7 +200,7 @@ def compute_inner_products(self, numerical=True, timeout=None, num_threads=None)
             eq.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
 
     def compute_tensor(self, numerical=True, compute_inner_products=False, compute_inner_products_kwargs=None,
-                       substitutions=None, basis_subs=False, parameters_subs=None):
+                       substitutions=None, basis_subs=False, parameters_subs=None, lhs_inversion=True):
         """Compute the tensor of the symbolic or numerical representation of the ordinary differential
         equations tendencies of the layer.
         Results are stored in the :attr:`~Layer.tensor` attribute.
@@ -212,7 +212,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             Whether to compute the numerical or the symbolic tensor.
             Default to `True` (numerical tensor as output).
         compute_inner_products: bool, optional
-            Whether the inner products tensors of the layer equations' terms must be compute first.
+            Whether the inner products tensors of the layer equations' terms must be computed first.
             Default to `False`. Please note that if the inner products are not computed firsthand, the tensor computation
             will fail.
         compute_inner_products_kwargs: dict, optional
@@ -229,6 +229,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
         parameters_subs: list(~parameter.Parameter), optional
             List of model's parameters to substitute in the symbolic tendencies' tensor.
             Only applies for the symbolic tendencies.
+        lhs_inversion: bool, optional
+            Try to inverse the LHS matrix and take the matricial product with the RHS. Default to `True`.
 
         """
 
@@ -260,7 +262,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             for field, eq in zip(self.fields, self.equations):
                 ndim = field.state.__len__()
                 if eq.other_fields_in_lhs:
-                    if self.other_fields_in_lhs:
+                    if self.other_fields_in_lhs or not lhs_inversion:
                         self._lhs_inversion = False
                         if self._cake is None:
                             raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
@@ -376,7 +378,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                             b_subs.append(sbsb)
                 ndim = field.state.__len__()
                 if eq.other_fields_in_lhs:
-                    if self.other_fields_in_lhs:
+                    if self.other_fields_in_lhs or not lhs_inversion:
                         self._lhs_inversion = False
                         if self._cake is None:
                             raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')

From 70191355ead8da6db2320467fc4eb600473a36f6 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 11:21:14 +0100
Subject: [PATCH 21/26] Solving LHS inversion logic

---
 layercake/bakery/cake.py   |  20 +++----
 layercake/bakery/layers.py | 108 ++++++++++++++++++-------------------
 2 files changed, 63 insertions(+), 65 deletions(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index 4eb6da2..8166f83 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -52,7 +52,7 @@ class Cake(object):
     def __init__(self):
 
         self.layers = list()
-        self._lhs_inversion = True
+        self._lhs_inversion_in_layer = True
 
     def add_layer(self, layer):
         """Add a layer object to the cake.
@@ -158,10 +158,10 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             Only applies for the symbolic tendencies.
 
         """
-        self._lhs_inversion = True
+        self._lhs_inversion_in_layer = True
         for layer in self.layers:
-            if not layer._lhs_inversion:
-                self._lhs_inversion = False
+            if layer.other_fields_in_lhs:
+                self._lhs_inversion_in_layer = False
                 break
 
         for layer in self.layers:
@@ -171,7 +171,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                  substitutions=substitutions,
                                  basis_subs=basis_subs,
                                  parameters_subs=parameters_subs,
-                                 lhs_inversion=self._lhs_inversion
+                                 lhs_inversion=self._lhs_inversion_in_layer
                                  )
 
     @property
@@ -219,11 +219,11 @@ def tensor(self):
 
         if numerical:
             tensor = sp.zeros(shape, dtype=np.float64, format='dok')
-            if self._lhs_inversion:
+            if not self._lhs_inversion_in_layer:
                 lhs_mat = sp.zeros((self.ndim + 1, self.ndim + 1), dtype=np.float64, format='dok')
         else:
             tensor = MutableSparseNDimArray(iterable={}, shape=shape)
-            if self._lhs_inversion:
+            if not self._lhs_inversion_in_layer:
                 lhs_mat = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
 
         for i, layer in enumerate(self.layers):
@@ -243,17 +243,17 @@ def tensor(self):
             args = tuple(slices + zeros)
             if numerical:
                 tensor[args] = tensor[args] + layer.tensor.todense()[1:]
-                if self._lhs_inversion:
+                if not self._lhs_inversion_in_layer:
                     lhs_mat[args[:2]] = lhs_mat[args[:1]] + layer._lhs_mat.todense()[1:]
             else:
                 tensor[args] = tensor[args] + layer.tensor[1:]
-                if self._lhs_inversion:
+                if not self._lhs_inversion_in_layer:
                     lhs_mat[args[:2]] = lhs_mat[args[:2]] + layer._lhs_mat[1:]
 
         if numerical:
             if self._lhs_inversion:
                 try:
-                    tensor = sp.COO(np.tensordot(np.linalg.inv(lhs_mat), tensor.to_coo(), 1))
+                    tensor = sp.COO(np.tensordot(np.linalg.inv(lhs_mat.todense()), tensor.to_coo(), 1))
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the cake is not invertible with the provided basis.')
             else:
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 39aa00d..0d1185f 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -200,7 +200,7 @@ def compute_inner_products(self, numerical=True, timeout=None, num_threads=None)
             eq.compute_inner_products(field.basis, numerical=numerical, timeout=timeout, num_threads=num_threads)
 
     def compute_tensor(self, numerical=True, compute_inner_products=False, compute_inner_products_kwargs=None,
-                       substitutions=None, basis_subs=False, parameters_subs=None, lhs_inversion=True):
+                       substitutions=None, basis_subs=False, parameters_subs=None, lhs_inversion_in_layer=True):
         """Compute the tensor of the symbolic or numerical representation of the ordinary differential
         equations tendencies of the layer.
         Results are stored in the :attr:`~Layer.tensor` attribute.
@@ -229,7 +229,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
         parameters_subs: list(~parameter.Parameter), optional
             List of model's parameters to substitute in the symbolic tendencies' tensor.
             Only applies for the symbolic tendencies.
-        lhs_inversion: bool, optional
+        lhs_inversion_in_layer: bool, optional
             Try to inverse the LHS matrix and take the matricial product with the RHS. Default to `True`.
 
         """
@@ -261,30 +261,29 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
             lhs_order = 1
             for field, eq in zip(self.fields, self.equations):
                 ndim = field.state.__len__()
-                if eq.other_fields_in_lhs:
-                    if self.other_fields_in_lhs or not lhs_inversion:
-                        self._lhs_inversion = False
-                        if self._cake is None:
-                            raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
-                        else:
-                            if np.all(self._lhs_mat.todense() == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                                self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
-                            if np.any(self._lhs_mat.todense() != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                                raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
-
-                        for lhs_term in eq.lhs_terms:
-                            ofield = lhs_term.field
-                            ofield_order = 1
-                            ondim = ofield.state.__len__()
-                            for tfield in self._cake.fields:
-                                if ofield is tfield:
-                                    break
-                                tndim = tfield.state.__len__()
-                                ofield_order += tndim
-                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] = \
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
-
+                if self.other_fields_in_lhs or not lhs_inversion_in_layer:
+                    self._lhs_inversion = False
+                    if self._cake is None:
+                        raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
                     else:
+                        if np.all(self._lhs_mat.todense() == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                            self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
+                        if np.any(self._lhs_mat.todense() != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                            raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+
+                    for lhs_term in eq.lhs_terms:
+                        ofield = lhs_term.field
+                        ofield_order = 1
+                        ondim = ofield.state.__len__()
+                        for tfield in self._cake.fields:
+                            if ofield is tfield:
+                                break
+                            tndim = tfield.state.__len__()
+                            ofield_order += tndim
+                        self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] = \
+                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
+
+                elif eq.other_fields_in_lhs:
                         for lhs_term in eq.lhs_terms:
                             ofield = lhs_term.field
                             ofield_order = 1
@@ -377,38 +376,37 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                         else:
                             b_subs.append(sbsb)
                 ndim = field.state.__len__()
-                if eq.other_fields_in_lhs:
-                    if self.other_fields_in_lhs or not lhs_inversion:
-                        self._lhs_inversion = False
-                        if self._cake is None:
-                            raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
-                        else:
-                            if np.all(self._lhs_mat == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                                self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
-                            if np.any(self._lhs_mat != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                                raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
-
-                        for lhs_term in eq.lhs_terms:
-                            ofield = lhs_term.field
-                            ofield_order = 1
-                            ondim = ofield.state.__len__()
-                            for tfield in self._cake.fields:
-                                if ofield is tfield:
-                                    break
-                                tndim = tfield.state.__len__()
-                                ofield_order += tndim
-                                self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
+                if self.other_fields_in_lhs or not lhs_inversion_in_layer:
+                    self._lhs_inversion = False
+                    if self._cake is None:
+                        raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
                     else:
-                        for lhs_term in eq.lhs_terms:
-                            ofield = lhs_term.field
-                            ofield_order = 1
-                            ondim = ofield.state.__len__()
-                            for tfield in self.fields:
-                                if ofield is tfield:
-                                    break
-                                tndim = tfield.state.__len__()
-                                ofield_order += tndim
+                        if np.all(self._lhs_mat == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                            self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
+                        if np.any(self._lhs_mat != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                            raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+
+                    for lhs_term in eq.lhs_terms:
+                        ofield = lhs_term.field
+                        ofield_order = 1
+                        ondim = ofield.state.__len__()
+                        for tfield in self._cake.fields:
+                            if ofield is tfield:
+                                break
+                            tndim = tfield.state.__len__()
+                            ofield_order += tndim
                             self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
+                elif eq.other_fields_in_lhs:
+                    for lhs_term in eq.lhs_terms:
+                        ofield = lhs_term.field
+                        ofield_order = 1
+                        ondim = ofield.state.__len__()
+                        for tfield in self.fields:
+                            if ofield is tfield:
+                                break
+                            tndim = tfield.state.__len__()
+                            ofield_order += tndim
+                        self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
                 else:
                     try:
                         lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = eq.lhs_inner_products_addition.inv().simplify()

From 2bdd6d4b4430e27f154bfff2e46342535a69026d Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 11:22:14 +0100
Subject: [PATCH 22/26] Cont'd

---
 layercake/bakery/cake.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index 8166f83..7840202 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -171,7 +171,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                  substitutions=substitutions,
                                  basis_subs=basis_subs,
                                  parameters_subs=parameters_subs,
-                                 lhs_inversion=self._lhs_inversion_in_layer
+                                 lhs_inversion_in_layer=self._lhs_inversion_in_layer
                                  )
 
     @property

From ad0693ac73384014812b53213d15c19f424c5d5e Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 11:29:12 +0100
Subject: [PATCH 23/26] Solved a bug and added a parameter to force LHS
 inversion at the cake level

---
 layercake/bakery/cake.py | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index 7840202..b114457 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -128,7 +128,7 @@ def number_of_layers(self):
         return self.layers.__len__()
 
     def compute_tensor(self, numerical=True, compute_inner_products=False, compute_inner_products_kwargs=None,
-                       substitutions=None, basis_subs=False, parameters_subs=None):
+                       substitutions=None, basis_subs=False, parameters_subs=None, lhs_inversion_in_layer=True):
         """Compute the tensor of the symbolic or numerical representation of the ordinary differential
         equations tendencies of all the layers.
         Arguments are passed to the layer :meth:`~Layer.compute_tensor` method.
@@ -156,13 +156,19 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
         parameters_subs: list(~parameter.Parameter), optional
             List of model's parameters to substitute in the symbolic tendencies' tensor.
             Only applies for the symbolic tendencies.
+        lhs_inversion_in_layer: bool, optional
+            Try to inverse the LHS matrix and take the matricial product with the RHS in each layer if possible. Default to `True`.
+            If `False`, it forces the inversion of the LHS at the cake level.
 
         """
-        self._lhs_inversion_in_layer = True
-        for layer in self.layers:
-            if layer.other_fields_in_lhs:
-                self._lhs_inversion_in_layer = False
-                break
+        if lhs_inversion_in_layer:
+            self._lhs_inversion_in_layer = True
+            for layer in self.layers:
+                if layer.other_fields_in_lhs:
+                    self._lhs_inversion_in_layer = False
+                    break
+        else:
+            self._lhs_inversion_in_layer = False
 
         for layer in self.layers:
             layer.compute_tensor(numerical=numerical,
@@ -251,7 +257,7 @@ def tensor(self):
                     lhs_mat[args[:2]] = lhs_mat[args[:2]] + layer._lhs_mat[1:]
 
         if numerical:
-            if self._lhs_inversion:
+            if not self._lhs_inversion_in_layer:
                 try:
                     tensor = sp.COO(np.tensordot(np.linalg.inv(lhs_mat.todense()), tensor.to_coo(), 1))
                 except LinAlgError:
@@ -259,7 +265,7 @@ def tensor(self):
             else:
                 tensor = tensor.to_coo()
         else:
-            if self._lhs_inversion:
+            if not self._lhs_inversion_in_layer:
                 try:
                     tensor = ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat.inv().simplify(), tensor, 1))
                 except NonInvertibleMatrixError:

From d7ab33bd43a6bd95bcf406966f20a83918e109dc Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Tue, 10 Mar 2026 14:31:41 +0100
Subject: [PATCH 24/26] Fixed bugs in the cake LHS inversion

---
 layercake/bakery/cake.py   | 8 ++++++--
 layercake/bakery/layers.py | 2 +-
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index b114457..c3197e4 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -259,7 +259,9 @@ def tensor(self):
         if numerical:
             if not self._lhs_inversion_in_layer:
                 try:
-                    tensor = sp.COO(np.tensordot(np.linalg.inv(lhs_mat.todense()), tensor.to_coo(), 1))
+                    lhs_mat_inverted = np.zeros((self.ndim + 1, self.ndim + 1))
+                    lhs_mat_inverted[1:, 1:] = np.linalg.inv(lhs_mat.todense()[1:, 1:])
+                    tensor = sp.COO(np.tensordot(lhs_mat_inverted, tensor.to_coo(), 1))
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the cake is not invertible with the provided basis.')
             else:
@@ -267,7 +269,9 @@ def tensor(self):
         else:
             if not self._lhs_inversion_in_layer:
                 try:
-                    tensor = ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat.inv().simplify(), tensor, 1))
+                    lhs_mat_inverted = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self.ndim + 1))
+                    lhs_mat_inverted[1:, 1:] = lhs_mat[1:, 1:].inv().simplify()
+                    tensor = ImmutableSparseNDimArray(symbolic_tensordot(lhs_mat_inverted, tensor, 1))
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the cake is not invertible with the provided basis.')
             else:
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 0d1185f..e04eaf5 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -480,7 +480,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
 
             if self._lhs_inversion and not self._lhs_inverted:
                 try:
-                    lhs_mat_inverted = self._lhs_mat.inv().simplify()
+                    lhs_mat_inverted[1:, 1:] = self._lhs_mat[1:, 1:].inv().simplify()
                     self._lhs_inverted = True
                 except NonInvertibleMatrixError:
                     raise NonInvertibleMatrixError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')

From 26d2fa053a837fc7fb5c67e59b707c04db64404e Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Wed, 11 Mar 2026 09:52:56 +0100
Subject: [PATCH 25/26] Corrected a bug due to weird sympy matrix indexing

---
 layercake/bakery/cake.py   | 2 +-
 layercake/bakery/layers.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py
index c3197e4..7b9e4e4 100644
--- a/layercake/bakery/cake.py
+++ b/layercake/bakery/cake.py
@@ -254,7 +254,7 @@ def tensor(self):
             else:
                 tensor[args] = tensor[args] + layer.tensor[1:]
                 if not self._lhs_inversion_in_layer:
-                    lhs_mat[args[:2]] = lhs_mat[args[:2]] + layer._lhs_mat[1:]
+                    lhs_mat[args[:2]] = lhs_mat[args[:2]] + layer._lhs_mat[1:, :]
 
         if numerical:
             if not self._lhs_inversion_in_layer:
diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index e04eaf5..82b4142 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -395,7 +395,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                 break
                             tndim = tfield.state.__len__()
                             ofield_order += tndim
-                            self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
+                        self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
                 elif eq.other_fields_in_lhs:
                     for lhs_term in eq.lhs_terms:
                         ofield = lhs_term.field

From ad51c88ff9fce8e01ecb90622ba8118f13e70783 Mon Sep 17 00:00:00 2001
From: Jonathan Demaeyer <jodemaey@meteo.be>
Date: Wed, 11 Mar 2026 11:56:50 +0100
Subject: [PATCH 26/26] Dealing again with errors and logic

---
 layercake/bakery/layers.py | 20 +++++++++-----------
 1 file changed, 9 insertions(+), 11 deletions(-)

diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py
index 82b4142..90f2f6e 100644
--- a/layercake/bakery/layers.py
+++ b/layercake/bakery/layers.py
@@ -264,12 +264,9 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 if self.other_fields_in_lhs or not lhs_inversion_in_layer:
                     self._lhs_inversion = False
                     if self._cake is None:
-                        raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
-                    else:
-                        if np.all(self._lhs_mat.todense() == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                            self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
-                        if np.any(self._lhs_mat.todense() != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                            raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+                        raise ValueError('Field(s) in the LHS are not found in the layer and there is no cake.')
+                    if self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                        self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
 
                     for lhs_term in eq.lhs_terms:
                         ofield = lhs_term.field
@@ -280,6 +277,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                 break
                             tndim = tfield.state.__len__()
                             ofield_order += tndim
+                        else:
+                            raise ValueError(f'Field {ofield} not found in the cake.')
                         self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] = \
                             self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
 
@@ -380,11 +379,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                     self._lhs_inversion = False
                     if self._cake is None:
                         raise ValueError('Field(s) in the LHS are not found in the layer or in the cake.')
-                    else:
-                        if np.all(self._lhs_mat == 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                            self._lhs_mat = sp.zeros((self.ndim + 1, self._cake.ndim + 1), dtype=np.float64, format='dok')
-                        if np.any(self._lhs_mat != 0) and self._lhs_mat.shape[1] != self._cake.ndim + 1:
-                            raise ValueError(f'Error in the initialization of the LHS in layer {self}.')
+                    if self._lhs_mat.shape[1] != self._cake.ndim + 1:
+                        self._lhs_mat = MutableSparseMatrix(sympy_zeros(self.ndim + 1, self._cake.ndim + 1))
 
                     for lhs_term in eq.lhs_terms:
                         ofield = lhs_term.field
@@ -406,6 +402,8 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                 break
                             tndim = tfield.state.__len__()
                             ofield_order += tndim
+                        else:
+                            raise ValueError(f'Field {ofield} not found in the cake.')
                         self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] += lhs_term.inner_products
                 else:
                     try: