@@ -26,7 +26,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal
 | [Neural Network](doc/neural_network/README.md)                     | Layers, activations, losses, model                                   |
 | [Optimization](doc/optimization/README.md)                         | Gradient Descent                                                     |
 | [Regularization](doc/regularization/README.md)                     | L1 (Lasso), L2 (Ridge)                                              |
-| [Math](doc/math/README.md)                                         | CORDIC, Quaternion, MatrixNorms, Step Response Metrics, MatrixExponential                                |
+| [Math](doc/math/README.md)                                         | CORDIC, Quaternion, MatrixNorms, Step Response Metrics, MatrixExponential                                                   |
 | [Solvers](doc/solvers/README.md)                                   | Gaussian Elimination, Levinson-Durbin, Durand-Kerner, Cholesky, DARE, Runge-Kutta ODE Integrators (RK4 + Dormand-Prince), Spectral Radius & Discrete Stability Margin, QR Decomposition (Householder / Givens), LU Decomposition with Partial Pivoting |
 | [Performance Optimization](doc/performance-optimization/README.md) | Compiler optimizations, SIMD                                         |
 
 
@@ -28,7 +28,6 @@ Difficulty legend:
 | #  | Component                                            | Target module             | Difficulty |
 |----|------------------------------------------------------|---------------------------|------------|
 | 30 | Continuous → discrete conversion (`c2d`)             | `math`                    | ★★★★☆      |
-| 31 | Lyapunov / Sylvester equation solvers                | `solvers`                 | ★★★★☆      |
 | 32 | Transfer-function ↔ state-space conversion           | `control_analysis`        | ★★★★☆      |
 | 33 | Madgwick / Mahony AHRS                               | `filters/active`          | ★★★★☆      |
 | 34 | Sliding Mode Control (SMC)                           | `robust_control` (new)    | ★★★★☆      |
 
@@ -0,0 +1,139 @@
+# Lyapunov and Sylvester Equation Solvers
+
+## Overview & Motivation
+
+The Sylvester equation $AX + XB = C$ and its symmetric specialisations — the continuous Lyapunov
+$AX + XA^\top = -Q$ and the discrete Lyapunov $AXA^\top - X = -Q$ — are the central matrix
+equations of stability analysis, gramian computation, and robust-control synthesis. A
+positive-definite Lyapunov solution is a formal certificate that the system $\dot{x} = Ax$ is
+asymptotically stable; the same equations produce the **controllability and observability gramians**
+that measure how effectively a linear system can be driven and observed.
+
+## Mathematical Theory
+
+### Sylvester Equation
+
+Given matrices $A \in \mathbb{R}^{N \times N}$, $B \in \mathbb{R}^{M \times M}$, and
+$C \in \mathbb{R}^{N \times M}$, the Sylvester equation is
+
+$$AX + XB = C,$$
+
+with unknown $X \in \mathbb{R}^{N \times M}$.  The equation has a unique solution if and only if
+$A$ and $-B$ share no eigenvalue, equivalently $\lambda_i(A) + \lambda_j(B) \neq 0$ for all
+$i, j$.
+
+### Kronecker Vectorisation
+
+Stacking $X$ column-by-column into $\operatorname{vec}(X) \in \mathbb{R}^{NM}$ and applying
+the identities $\operatorname{vec}(AXI) = (I_M \otimes A)\operatorname{vec}(X)$ and
+$\operatorname{vec}(IXB) = (B^\top \otimes I_N)\operatorname{vec}(X)$ converts the matrix
+equation into the standard linear system
+
+$$\underbrace{(I_M \otimes A + B^\top \otimes I_N)}_{K \in \mathbb{R}^{NM \times NM}}
+\operatorname{vec}(X) = \operatorname{vec}(C).$$
+
+The $(iN+k, jN+l)$ entry of $K$ is $\delta_{ij} A_{kl} + B_{ji} \delta_{kl}$.  The system
+is solved by Gaussian elimination with partial pivoting; a near-zero pivot signals the
+unsolvable (degenerate eigenvalue) case.
+
+### Continuous Lyapunov Equation
+
+Setting $B = A^\top$ and $C = -Q$ reduces the Sylvester equation to the continuous Lyapunov form
+
+$$AX + XA^\top = -Q.$$
+
+When $A$ is Hurwitz (all eigenvalues in the open left half-plane) and $Q$ is symmetric
+positive semi-definite, the unique solution $X$ is symmetric positive definite.
+
+### Discrete Lyapunov Equation
+
+The discrete counterpart
+
+$$AXA^\top - X = -Q$$
+
+is solved by the Kronecker system $(A \otimes A - I_{N^2})\operatorname{vec}(X) = -\operatorname{vec}(Q)$,
+with the $(iN+k, jN+l)$ entry of $A \otimes A - I$ equal to $A_{ij}A_{kl} - \delta_{ij}\delta_{kl}$.
+When $A$ is Schur-stable ($\rho(A) < 1$) and $Q \succeq 0$, the unique solution equals the
+controllability gramian $X = \sum_{k=0}^{\infty} A^k Q (A^\top)^k$.
+
+## Complexity Analysis
+
+| Case                          | Time           | Space                       | Notes                                                         |
+|-------------------------------|----------------|-----------------------------|---------------------------------------------------------------|
+| $N{=}M{=}2$                   | $O(64)$ flops  | $O(NM)^2 = 16$ floats for K | 4×4 linear system                                             |
+| $N{=}M{=}3$                   | $O(729)$ flops | 81 floats for K             | 9×9 linear system                                             |
+| General $N{=}M$               | $O(N^6)$ flops | $O(N^4)$ words              | Gaussian elimination on the $N^2 \times N^2$ Kronecker system |
+| Alternative (Bartels–Stewart) | $O(N^3)$       | $O(N^2)$                    | Requires Schur decomposition; preferred for $N \geq 6$        |
+
+For the embedded sizes ($N \leq 5$) targeted by this library the Kronecker approach is negligibly
+more expensive than Bartels–Stewart and avoids the full Schur machinery.
+
+## Step-by-Step Walkthrough
+
+Solve $AX + XB = C$ for the $2 \times 2$ case
+$A = \bigl[\begin{smallmatrix}-3&1\\0&-2\end{smallmatrix}\bigr]$,
+$B = \bigl[\begin{smallmatrix}1&2\\0&1\end{smallmatrix}\bigr]$,
+$C = I_2$.
+
+**Build $K$** (the $4 \times 4$ Kronecker matrix):
+
+$$K = \begin{pmatrix}-2&1&0&0\\0&-1&0&0\\2&0&-2&1\\0&2&0&-1\end{pmatrix}.$$
+
+**Build $\operatorname{vec}(C) = (1, 0, 0, 1)^\top$.**
+
+**Gaussian elimination** with partial pivoting produces $X_{00} = -0.5$, $X_{10} = 0$,
+$X_{01} = -1$, $X_{11} = -1$.
+
+**Verify:** $AX + XB = I_2$. ✓
+
+## Pitfalls & Edge Cases
+
+- **Unsolvable case:** if $\lambda_i(A) + \lambda_j(B) = 0$ (continuous) or
+  $\lambda_i(A)\lambda_j(B) = 1$ (discrete) for any pair, $K$ is singular and no unique solution
+  exists.  The solver returns `false` and leaves the solution unchanged.
+- **Ill-conditioning:** near-degenerate eigenvalue pairs yield large solution norms and reduced
+  accuracy.  Monitor the residual $\|AX + XB - C\|$ after solving.
+- **Symmetry:** the continuous Lyapunov solution is guaranteed symmetric only when $Q$ is
+  symmetric; numerical symmetry should be enforced by the caller if required downstream.
+
+## Variants & Generalizations
+
+- **Bartels–Stewart algorithm:** reduces to quasi-upper-triangular (Schur) form via the double-shift
+  QR algorithm, then back-substitutes column by column in $O(N^3)$.  Preferred for $N \geq 6$.
+- **Smith iteration:** fixed-point iteration $X_{k+1} = AX_k A^\top + Q$ for the discrete
+  Lyapunov equation; simple but slow for $\rho(A)$ close to 1.
+- **Cholesky factor updates:** propagate a Cholesky factor of $X$ for guaranteed positive-definiteness
+  in square-root filtering (see `SquareRootKalmanFilter`).
+
+## Applications
+
+- **LQR / LQG design:** the continuous Lyapunov equation $AX + XA^\top + Q = 0$ yields the
+  closed-loop covariance, and the DARE solution can be cross-validated via its Lyapunov residual.
+- **Kalman filter (observability Gramian):** $AW + WA^\top + BB^\top = 0$ gives the
+  controllability Gramian used to evaluate reachability in balanced truncation.
+- **Pole-placement observers (Luenberger):** the Sylvester equation $AX - XF = LC$ arises when
+  placing observer poles while decoupling the error dynamics.
+- **H₂ norm computation:** $\|G\|_2^2 = \operatorname{tr}(B^\top X B)$ where $X$ satisfies the
+  observability Gramian Lyapunov equation.
+
+## Connections to Other Algorithms
+
+- **LU Decomposition:** the Kronecker-vectorisation approach converts any Sylvester equation into
+  a dense linear system solved by the `LuDecomposition` solver.
+- **QR Decomposition:** the Bartels–Stewart and Hammarling algorithms use QR / Schur
+  decomposition to reduce the problem to quasi-triangular form before back-substitution.
+- **DARE:** the discrete algebraic Riccati equation can be warm-started or verified using
+  one step of the discrete Lyapunov equation.
+- **Spectral Radius:** solvability of the continuous equation requires $\lambda_i(A) + \lambda_j(A) \neq 0$;
+  checking the spectral radius is a fast necessary condition.
+
+## References & Further Reading
+
+- Bartels, R. H. & Stewart, G. W. (1972). "Solution of the matrix equation AX + XB = C."
+  *Communications of the ACM*, 15(9), 820–826.
+- Golub, G. H., Nash, S. & Van Loan, C. (1979). "A Hessenberg-Schur method for the problem
+  AX + XB = C." *IEEE Transactions on Automatic Control*, 24(6), 909–913.
+- Horn, R. A. & Johnson, C. R. (1994). *Topics in Matrix Analysis*. Cambridge University Press.
+  (Chapter 4: Kronecker products and the vec operator.)
+- Hammarling, S. J. (1982). "Numerical solution of the stable, non-negative definite Lyapunov
+  equation." *IMA Journal of Numerical Analysis*, 2(3), 303–323.
@@ -16,3 +16,4 @@ Numerical solvers for linear systems, polynomial roots, and matrix equations.
 | [Condition Number](ConditionNumber.md)                                     | 1-norm condition number estimate via column-wise inverse solve         |
 | [QR Decomposition](QrDecomposition.md)                                     | Householder factorization and Givens streaming row update for least-squares solves |
 | [LU Decomposition](LuDecomposition.md)                                     | PA = LU factorization with partial pivoting for general dense linear systems       |
+| [Lyapunov / Sylvester Solvers](LyapunovSylvester.md)                       | Sylvester AX+XB=C and continuous/discrete Lyapunov solvers via Kronecker vectorisation |
@@ -19,6 +19,7 @@ target_sources(numerical.solver PRIVATE
     GaussianElimination.hpp
     LevinsonDurbin.hpp
     LuDecomposition.hpp
+    LyapunovSylvester.hpp
     OdeSystem.hpp
     QrDecomposition.hpp
     RungeKuttaIntegrators.hpp
@@ -31,6 +32,7 @@ numerical_add_coverage_sources(numerical.solver
     DurandKerner.cpp
     GaussianElimination.cpp
     LuDecomposition.cpp
+    LyapunovSylvester.cpp
     QrDecomposition.cpp
     RungeKuttaIntegrators.cpp
     SpectralRadius.cpp
 
@@ -0,0 +1,7 @@
+#include "numerical/solvers/LyapunovSylvester.hpp"
+
+namespace solvers
+{
+    template class SylvesterSolver<float, 2, 2>;
+    template class SylvesterSolver<float, 3, 3>;
+}
@@ -0,0 +1,148 @@
+#pragma once
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC optimize("O3", "fast-math")
+#endif
+
+#include "numerical/math/CompilerOptimizations.hpp"
+#include "numerical/math/Matrix.hpp"
+#include "numerical/solvers/LuDecomposition.hpp"
+#include <cstddef>
+#include <type_traits>
+
+namespace solvers
+{
+    template<typename T, std::size_t N, std::size_t M>
+    class SylvesterSolver
+    {
+        static_assert(std::is_floating_point_v<T>, "SylvesterSolver supports floating-point types only");
+        static_assert(N > 0 && M > 0, "SylvesterSolver requires positive dimensions");
+
+        static constexpr std::size_t NM = N * M;
+
+    public:
+        SylvesterSolver() = default;
+
+        OPTIMIZE_FOR_SPEED bool SolveSylvester(
+            const math::SquareMatrix<T, N>& a,
+            const math::SquareMatrix<T, M>& b,
+            const math::Matrix<T, N, M>& c);
+
+        OPTIMIZE_FOR_SPEED bool SolveContinuousLyapunov(
+            const math::SquareMatrix<T, N>& a,
+            const math::SquareMatrix<T, N>& q);
+
+        OPTIMIZE_FOR_SPEED bool SolveDiscreteLyapunov(
+            const math::SquareMatrix<T, N>& a,
+            const math::SquareMatrix<T, N>& q);
+
+        [[nodiscard]] const math::Matrix<T, N, M>& Solution() const
+        {
+            return solution;
+        }
+
+    private:
+        static math::Matrix<T, N * M, 1> Vec(const math::Matrix<T, N, M>& x);
+        static math::Matrix<T, N, M> UnVec(const math::Matrix<T, N * M, 1>& v);
+
+        math::Matrix<T, N, M> solution{};
+    };
+
+    template<typename T, std::size_t N, std::size_t M>
+    OPTIMIZE_FOR_SPEED bool SylvesterSolver<T, N, M>::SolveSylvester(
+        const math::SquareMatrix<T, N>& a,
+        const math::SquareMatrix<T, M>& b,
+        const math::Matrix<T, N, M>& c)
+    {
+        math::Matrix<T, NM, NM> K{};
+        for (std::size_t row = 0; row < NM; ++row)
+            for (std::size_t col = 0; col < NM; ++col)
+            {
+                const std::size_t i = row / N;
+                const std::size_t k = row % N;
+                const std::size_t j = col / N;
+                const std::size_t l = col % N;
+                T val{};
+                if (i == j)
+                    val += a.at(k, l);
+                if (k == l)
+                    val += b.at(j, i);
+                K.at(row, col) = val;
+            }
+
+        math::Matrix<T, NM, 1> vecC = Vec(c);
+
+        LuDecomposition<T, NM> lu{};
+        if (!lu.Decompose(K))
+            return false;
+
+        solution = UnVec(lu.Solve(vecC));
+        return true;
+    }
+
+    template<typename T, std::size_t N, std::size_t M>
+    OPTIMIZE_FOR_SPEED bool SylvesterSolver<T, N, M>::SolveContinuousLyapunov(
+        const math::SquareMatrix<T, N>& a,
+        const math::SquareMatrix<T, N>& q)
+    {
+        static_assert(N == M, "SolveContinuousLyapunov requires square solver (N == M)");
+
+        return SolveSylvester(a, a.Transpose(), q * T(-1));
+    }
+
+    template<typename T, std::size_t N, std::size_t M>
+    OPTIMIZE_FOR_SPEED bool SylvesterSolver<T, N, M>::SolveDiscreteLyapunov(
+        const math::SquareMatrix<T, N>& a,
+        const math::SquareMatrix<T, N>& q)
+    {
+        static_assert(N == M, "SolveDiscreteLyapunov requires square solver (N == M)");
+
+        math::Matrix<T, NM, NM> K{};
+        for (std::size_t row = 0; row < NM; ++row)
+            for (std::size_t col = 0; col < NM; ++col)
+            {
+                const std::size_t i = row / N;
+                const std::size_t k = row % N;
+                const std::size_t j = col / N;
+                const std::size_t l = col % N;
+                T val = a.at(i, j) * a.at(k, l);
+                if (i == j && k == l)
+                    val -= T(1);
+                K.at(row, col) = val;
+            }
+
+        math::Matrix<T, NM, 1> rhs = Vec(q * T(-1));
+
+        LuDecomposition<T, NM> lu{};
+        if (!lu.Decompose(K))
+            return false;
+
+        solution = UnVec(lu.Solve(rhs));
+        return true;
+    }
+
+    template<typename T, std::size_t N, std::size_t M>
+    math::Matrix<T, N * M, 1> SylvesterSolver<T, N, M>::Vec(const math::Matrix<T, N, M>& x)
+    {
+        math::Matrix<T, NM, 1> v{};
+        for (std::size_t i = 0; i < M; ++i)
+            for (std::size_t k = 0; k < N; ++k)
+                v.at(i * N + k, 0) = x.at(k, i);
+        return v;
+    }
+
+    template<typename T, std::size_t N, std::size_t M>
+    math::Matrix<T, N, M> SylvesterSolver<T, N, M>::UnVec(const math::Matrix<T, N * M, 1>& v)
+    {
+        math::Matrix<T, N, M> out{};
+        for (std::size_t i = 0; i < M; ++i)
+            for (std::size_t k = 0; k < N; ++k)
+                out.at(k, i) = v.at(i * N + k, 0);
+        return out;
+    }
+
+#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
+    extern template class SylvesterSolver<float, 2, 2>;
+    extern template class SylvesterSolver<float, 3, 3>;
+#endif
+}
@@ -14,6 +14,7 @@ target_sources(numerical.solvers_test PRIVATE
     TestGaussianElimination.cpp
     TestLevinsonDurbin.cpp
     TestLuDecomposition.cpp
+    TestLyapunovSylvester.cpp
     TestQrDecomposition.cpp
     TestRungeKuttaIntegrators.cpp
     TestSpectralRadius.cpp