I have the following system of equations
$$ \sum_{l=1}^n D_{il}H_{lj}^i=\sum_{l=1}^n\sum_{m=1}^nR_{lm}^jD_{lm}+K_j $$ for all $i\in\{1,\dots,n\}$ and $j\in\{1,\dots,n\}$. I know the tensors $H$ and $R$, and the vector $K$. I would like to solve for $D$. Any help would be much appreciated!
Edit: I had forgotten the terms $K_j$.
This system is linear for $D_{i,j}$ variables. there are $n^2$ variables and $n^2$ equations so you can solve it with changing this to linear system. you should just rearrange matrix notations in $M_{n \times n}(\mathbb{R})$ to $\mathbb{R}^{n^2}$.
first make equation more suitable for factorization : $$ \sum_l H^i_{l,j} D_{il} - \sum_l \sum_m R^j_{l,m} D_{l,m} = K_j $$ raplace dummy index for first term: $$ \sum_m H^i_{m,j} D_{im} - \sum_l \sum_m R^j_{l,m} D_{l,m} = K_j $$ make a dummy summation over $l$ in first term: $$ \sum_l \sum_m H^i_{m,j} \delta_{i,l} D_{im} - \sum_l \sum_m R^j_{l,m} D_{l,m} = K_j $$ $$ \sum_l \sum_m \big( H^i_{m,j} \delta_{i,l} - R^j_{l,m} \big) D_{l,m} = K_j \tag{1} \label{1} $$ in which $\delta_{i,l}$ is Kronecker delta.
Now introduce new vector $\hat{D}$ for replacing matrix $D_{l,m}$. $$ \hat{D}_k := D_{l(k),m(k)} $$ where $l(k) = \lfloor \frac{k}{n} \rfloor +1$ and $m(k) = k-\lfloor \frac{k}{n} \rfloor n$. you can arrange differently, i just put every row of $D_{l,m}$ consecutive after each others. $\ref{1}$ becomes : $$ \sum_k \big( H^i_{m(k),j} \delta_{i,l(k)} - R^j_{l(k),m(k)} \big) \hat{D}_{k} = K_j .\tag{2}\label{2}$$ Now we should rearrange $(i,j)$ index to $1$ index for creating matrix system for $\ref{2}$. like before let $i(s) = \lfloor \frac{s}{n} \rfloor +1$ and $j(s) = s-\lfloor \frac{s}{n} \rfloor n$. so we obtain : $$ A \hat{D} = B \tag{3}\label{3} $$ with $A \in M_{n^2 \times n^2}(\mathbb{R})$ , $\hat{D},B \in \mathbb{R}^{n^2}$, in which : $$ A_{\alpha,\beta} = H^{i(\alpha)}_{m(\beta),j(\alpha)} \delta_{i(\alpha),l(\beta)} - R^{j(\alpha)}_{l(\beta),m(\beta)} \quad \alpha,\beta \in \{1,\dots,n^2 \} \\ B_{\alpha} = K_{j(\alpha)}. \tag{4}\label{4} $$ So if $\det(A) \neq 0$ the solution would be $\hat{D} = A^{-1} B$ then you can recover $D$ from $\hat{D}$.
Even for small $n$, $A$ is a large matrix so you should try solving this via a numerical computation software.
functions for $l(k)$ and $i(k)$ are same and also for $m(k)$ and $j(k)$ but i use this different notation for $\ref{4}$ not to be confusing. Also because of $\delta$ function you can also consider upper index of $H$ as $l$ from equation $\ref{1}$, it doesn't matter.