How to Solve a linear matrix equation of an array $M = BMC$ where $ B$ and $C$ are known Adding to the question's description : 
I am doing Feature extraction from videos and i am trying to implement this one line of mathematical equation to matlab or even any algorithm .
let's say I have $ B $ and $ C $ which are $ 10 ×10$ matrices 


*

*$ R $ is also a $ 10 × 10 $ matrix but it should satisfy two conditions* 


1) $M = B×M×C$
2) The Values in $ M $ should be within a given range [$x$( which is min) $y$(max)] given by the user 
i can handle the constraints of the range but i have no idea on how to proceed with the 1st condition at all , much less converting it into matlab
I have found several solutions to this type of problem
1)Sylvester's equation which approximates $ AX+XB = C $ or $AX + XB + C = 0 $
2)Lyapunov Matrix equation which approximates $AX+XA'+Q = 0$
both of these implementations of matlab i have found in this link given below 
https://github.com/ajt60gaibb/freeLYAP
but i need help in the mathematics area as to how i can use those (if i can) to  $M = BMC$
 A: If you multiply the equation $M = BMC$ by the inverse of $B$ on both sides and take everything on the same side of the equal sign, you get $$ B^{-1} M - MC = 0.$$
This is a Sylvester equation where the last term is the null matrix (all elements equal to zero). You should be able to solve your problem in Matlab by using lyap(inv(B),-C,zeros(10,10)) (here it should be ok to use the inv command because you are dealing with 10 x 10 matrices, consider another way of inverting $B$ if you also need to solve the same problem for bigger ones).
A: This equation is a special case of the discrete Sylvester equation, which takes the form $$X + AXB = C$$ for unknown matrix X. Solution $X$ is given by $$(I + B^⊤ \otimes A)vec(X) = vec(C)$$
where $vec(Q)$ is vectorization operator, that for a matrix $Q$ of size $m \times n$ creates a vector of size $mn\times 1$ by concatenating columns of Q, and $\otimes$ denotes the Kronecker product. This representation of the solution $X$ should not be used to solve this equation in general (there are much more efficient methods), but is useful to construct solution of the original problem.
In the original problem the matrix C is zero and solution X is given by
$$vec(X) = N x\\ N = null(I + B^⊤ \otimes A)$$
where $null(Q)$ is the null space of the matrix Q and $x$ is any vector of appropriate length. Therefore one must find any vector $x$ such that
$$a \leq Nx \leq b \;\;\;\;\; (LP)$$ where $a, b$ are lower and upper bounds on elements of the solution $X$.
In order to find such $x$ one must find any feasible solution to linear programming problem (LP).
