  A numerical method is developed leading to algebraic systems based on
generalized Lyapunov-Sylvester operators to approximate the solution of
two-dimensional Kuramoto-Sivashinsky equation. It consists of an order
reduction method and a finite difference discretization which is proved to be
uniquely solvable, stable and convergent by using Lyapunov criterion and
manipulating generalized Lyapunov-Sylvester operators. Some numerical
implementations are provided at the end to validate the theoretical results.
