Timo Hanke has written on this problem in the generality of cyclic algebras (http://arxiv.org/abs/math/0702681). He shows that it is equivalent to the solution of a norm equation, which has an algorithmic solution over global fields. In your case, this would in general involve the solution of a norm equation over a biquadratic field; I don't know how well this would fare in practice.
In the special case of quaternion algebras that you are inquiring about, there is a link to quadratic forms, and as Keith Conrad mentions, it is equivalent to find a zero of a quadratic form in six variables. This goes back to Albert, who looked at this form in detail in order to prove that there was an "honest" (non-quaternion) biquaternion algebra. A good reference for this is section 16 of the "Book of Involutions" or section XII.2 of Lam's "Introduction to Quadratic Forms over Fields". It may seem like just a reformulation of the problem, but it turns out there are algorithmic methods to find points on quadrics over number fields that are quite efficient in practice: the buzzword here is "indefinite LLL", and Watkins (http://magma.maths.usyd.edu.au/~watkins/papers/illl.pdf) explains what Magma does to accomplish this task.
In general, here is an idea I kicked around once which at least reduces the problem to solve norm equations over quadratic extensions (instead of biquadratic extensions). This might be only of theoretic/algorithmic interest, but at least it potentially generalizes. We wish to test if $A \cong B$ over a global field $F$ (say of characteristic not $2$ for now), and if so, to find an explicit isomorphism. If $A=(a,b)$ and $B=(c,d)$ and $a=c$, then there is an isomorphism if and only if $b/d$ is a norm from $\mathbb{Q}(\sqrt{a})$, and this can be accomplished algorithmically by a norm equation over this field; so it is enough to reduce to this case. To find a common subfield $K=\mathbb{Q}(\sqrt{a})$ in $A,B$, one can simply pick one (choose $K$ such that $K_v$ is not split at all places $v$ ramified in $A$ and $B$, e.g., take $a=-\mathrm{lcm}(ab,cd)$ if $\gcd(a,b)=\gcd(c,d)=1$, so this step does not even require factoring). The problem of embedding a quadratic field $K$ in a quaternion algebra is equivalent to (checking if $K$ splits the algebra and so) to a norm equation (a standard result, see e.g. http://www.math.dartmouth.edu/~jvoight/articles/quatalgs-060513.pdf). Once the field is embedded, we can diagonalize the quadratic form to reduce to the case where $a=c$.
This approach probably generalizes to cyclic algebras (of any characteristic), and if it is interesting to you, it is something I would be happy to work out with you.
