How to compute the pure extensions of a given mixed state?
The fact that every mixed state $\rho$ acting on a finite dimensional Hilbert spaces can be viewed as the reduced state of some pure state $|\psi\rangle$ on a bigger Hilbert space is known as purification, see this Wikipedia page, where also the algorithm is given.
In OP's case of $$\rho~\in~\mathcal{B}(\mathbb{C}^n\otimes \mathbb{C}^n),$$ one may choose a pure state $|\psi\rangle$ in the following Hilbert space $$|\psi\rangle~\in~\mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n.$$
It is a general result that any mixed quantum state can be viewed as the reduced state of a pure state in a larger dimensional Hilbert space. This is referred to as a purification, and some people even refer to the power of this idea as the "Church of the Larger Hilbert Space".
There exists a canonical way of constructing a purification, which has the advantage that it immediately demonstrates what the minimum dimension of the larger Hilbert space must be. So here is how you construct a purification: Let $\rho$ be a state in a Hilbert space of dimension $d$ with eigenvectors $\{|\phi_i\rangle\}_{i=1}^d$. Then
$|\psi\rangle=\sum_{i=1}^d|i\rangle\otimes|\phi_i\rangle$
is a purification of $\rho$, where $\{|i\rangle\}$ is an orthonormal set of vectors.
Therefore, you can see that it is always possible to construct a purification and moreover, the enlarged Hilbert space in general must have dimension at least the twice as big as the original one.
You correctly say that such a state won't be unique. Indeed, for example, pick the 50%-50% statistical mixture of two random pure basis states in ${\mathbb C}^n\otimes {\mathbb C}^n$. Each of these two states may be entangled with one of two orthogonal states in the third ${\mathbb C}^n$; however, what these two states in the third factor Hilbert space happen to be is completely undetermined.
So even if you had a constructive method to find the pure state in the three-part Hilbert space, it would fail to generate unique results.
However, there is another, much more serious problem with your proposal: in almost all cases, it has no solutions at all. Indeed, it's easy to demonstrate this fact by a simple counting of degrees of freedom. Pure states in ${\mathbb C}^n\otimes {\mathbb C}^n$ are specified by $n^2$ different complex numbers (one of them, the overall complex normalization, is unphysical).
Similarly, a general density matrix on this space is a $n^2\times n^2$ Hermitian matrix so it contains $n^4$ independent real parameters (one of them is the trace which should probably be set to one). However, that's larger than $2n^3$ number of real parameters coming from $n^3$ complex parameters of a wave function in ${\mathbb C}^n\otimes {\mathbb C}^n\otimes{\mathbb C}^n.$ At least for $n\gt 2$, it's larger. So up to a measure-zero subset of cases, you won't be able to find any pure state that reduces to the given mixed state. The diversity of the required results (density matrices) is much larger than the diversity of the ingredients (pure three-block states) that you may use to produce the desired outcome.
Of course, if you only had a density matrix for one of the three blocks, and not two of them, you would be able to solve it. At least the counting of the parameters wouldn't make the existence of a solution for a generic density matrix impossible.