I have been learning about the concept of entanglement swapping and found an equation mentioned in the textbook, Mathematics Of Quantum Computing: An Introduction written by Wolfgang Scherer.
At section 4.13, equation (4.13), $$$$\begin{split} {|\Phi\rangle}^{ABCD} & = |\Psi^{-}\rangle^{AB} \otimes |\Psi^{-}\rangle^{CD} \\ & = \frac{1}{2} \left( |0101\rangle - |0110\rangle -|1001\rangle + |1010\rangle \right)\\ & = \frac{1}{2} \left( |\Psi^{+}\rangle^{AD} \otimes |\Psi^{+}\rangle^{BC} - |\Psi^{-}\rangle^{AD} \otimes |\Psi^{-}\rangle^{BC}\\ \quad - |\Phi^{+}\rangle^{AD} \otimes |\Phi^{+}\rangle^{BC} + |\Phi^{-}\rangle^{AD} \otimes |\Phi^{-}\rangle^{BC} \right)\\ \end{split}$$$$ where
$$$$\begin{split} |\Phi^{+}\rangle = \frac{1}{\sqrt{2}} \left( |00\rangle + |11\rangle \right)\\ |\Phi^{-}\rangle = \frac{1}{\sqrt{2}} \left( |00\rangle - |11\rangle \right)\\ |\Psi^{+}\rangle = \frac{1}{\sqrt{2}} \left( |01\rangle + |10\rangle \right)\\ |\Psi^{-}\rangle = \frac{1}{\sqrt{2}} \left( |01\rangle - |10\rangle \right) \end{split}$$$$ and the supernotes $$\{ A,B,C,D \}$$ outside the ket indicate the qubits forming the bell state. Therefore, by doing bell measurement on qubits $$(B,C)$$, qubits $$(A,D)$$ are in bell state. This result is also mentioned in Multiparticle generalization of entanglement swapping by S. Bose et al (PhysRevA.57.822), but without derivation.

So my question is that how to derive other similar equations efficiently instead of expand the bell states back into computational basis $$| x \rangle^{\otimes n}$$ where $$x \in \{0,1\}$$ and go through the tedious algebras step by step? For example, how to do the entanglement swapping of $${|\Phi\rangle}^{ABCD} = |\Psi^{-}\rangle^{AB} \otimes |\Phi^{+}\rangle^{CD}$$ from $$(A,B)$$ and $$(C,D)$$ to $$(A,D)$$ and $$(B,C)$$ or even $$(A,C)$$ and $$(B,D)$$?

I guess the way that I'd start (aside from just getting a computer to do it!) is to remember that the Bell states form an orthonormal basis. So, you can ask, for example, about what the $$|\Phi^+\rangle^{AD}$$ component is: $$\langle\Phi^+|^{AD}|\Phi\rangle^{ABCD}=-\frac12|\Phi^+\rangle^{BC}.$$ You do this for each of the four states, and you can use that to reconstruct the overall description that you've given.
Now, for simplifying the above calculation without expanding out the state in the computational basis. Remember that all Bell states can be converted into $$|\Phi^+\rangle$$ by Pauli operations on either qubit. So, reduce the calculation to $$\langle\Phi^+|^{AD}|\Psi^-\rangle^{AB}|\Psi^-\rangle^{CD}=Y_BY_C\langle\Phi^+|^{AD}|\Phi^+\rangle^{AB}|\Phi^+\rangle^{CD}.$$ How does this help? Well, you could just explicitly calculate this once, or believe me, that $$\langle\Phi^+|^{AD}|\Phi^+\rangle^{AB}|\Phi^+\rangle^{CD}=\frac12|\Phi^+\rangle^{BC},$$ (the idea being that you can use this one result again and again) so the outcome is $$\frac12 Y_BY_C|\Phi^+\rangle^{BC}=-\frac12|\Phi^+\rangle^{BC}.$$