Solving equations involving both matrix and three-index unknowns Suppose I have two equations where the two unknowns and constants are square matrices. That's easy to solve since I can invert the matrices. But what if I have something like
$$A_{ik}X_{kj}+B_{ikl}Y_{klj}=\mathbb{I}_{ij}$$
$$C_{ijl}X_{lk}+D_{ijlm}Y_{lmk}=0_{ijk}$$
where one of the unknowns ($X$) is a matrix while the other ($Y$) is a three-index object (not necessarily a tensor) and the constants are two, three or four index objects properly contracted. Should I just give up at this point, or is there some way to go about "solving" this and expressing $X$ and $Y$ in terms of $A$, $B$, $C$ and $D$?
EDIT: I'm using summation notation (repeated indices in each term are summed over).
All the indices run from 1 to N, so you either have square matrices or NxNxN objects.
An algebraic solution is what I'm ideally looking for. As a start, I can express the first equation as $$X=A^{-1}(\mathbb{I}-BY)$$ since A is an (invertible) matrix. After substituting this into the second equation, I run into trouble solving for $Y$ because $D$ doesn't have a similar inversion I can utilize.
 A: If the objects in your equations are unstructured arrays then your equations are just a very concise way of expressing a potentially huge (depending on the range of indexes) system of linear equations on scalars. While this might be a neat way to write it, I do not see a whole lot you can do with this representation algebraically, unless there are some constraints or symmetries on the objects involved. If you have concrete values for the components of A, B, C & D use a computer to solve these systems numerically by generating the linear equations and using a linear solver.
However it looks like you will end up with many more variables than equations.
Also it looks like you are using the Einstein summation convention in these equations, since otherwise the equations would not be well formed.
This is not obvious at first though, since you are using only lower indexes and you say the entities are not tensors. You should probably explicitly add the summation signs or explicitly say you are using the Einstein convention and maybe raise one of each of the index pairs being summed. While the Einstein convention very useful in tensor calculus, I would lean towards explicit summation signs in your examples
A: Because you said the case of two unknown matrices is "easy to solve", I assume that you do not want to solve this by hand, but merely map it to a problem of a more usual form, allowing you to use numerical standard routines. 
You can reexpress your system by formally defining $Z_{ijk}:=X_{ik}$ for $j=1...N$ and it becomes
$$A_{ik}Z_{krj}+B_{ikl}Y_{klj}=\mathbb{I}_{ij}\quad (\forall ij)$$
$$C_{ijl}Z_{lrk}+D_{ijlm}Y_{lmk}=0_{ijk}\quad(\forall ijk)$$
Now, by definition the value of $Z_{ijk}$ must be independ of its second index, which we can assure by adding the equations
$$Z_{ijk}=Z_{i,j+1,k}\quad (\forall ijk,j\neq N)$$
Note that they are also linear, and including them we again arrive at the same number of equations and unknowns. Finally I did not tell you how to choose $N$, which is indeed free to you. You might for example want the second indices of $Y$ and $Z$ to be in the same range (which is already the case for the third index, as follows from the equations).
