One can solve the problem explicitly over $\mathbb{C}$ then try to work out the extra constraints due to working over the integers.
For the complex case, and continuing on the approach in the comment by Frank Thorne,
you can also use the so-called canonizant. Let $C(x,y)$ be your cubic.
The canonizant here is the same as the Hessian which classically normalized is:
$$
H(x,y)=\frac{1}{36}
\left(\frac{\partial^2}{\partial x \partial v}-\frac{\partial^2}{\partial y\partial u}\right)^2\ C(x,y) C(u,v)\ |_{u:=x, v:=y}
$$
meaning: take the derivatives then set $u=x$ and $v=y$.
If one can write $C=L_1^3+L_2^3$ where $L_1$, $L_2$ are linear forms in $x,y$,
then
$$
H(x,y)=2 \Delta(L_1,L_2)^2\ L_1(x,y)\ L_2(x,y)
$$
where $\Delta(L_1,L_2)$ is the determinant formed with the coefficients of the two linear forms.
The matrix you are looking for essentially is that which sends $(x,y)$
to $(L_1(x,y),L_2(x,y))$. So to find it you need to compute the Hessian and then factor it.
This means solving a second degree equation instead of a cubic equation as suggested in Frank's comment.