(* EPR = ((|01>-|10>)/(√2)) *)
EPR=(1/√2)*(({{0}, {1}, {-1}, {0}}))

(* Ψ = ((|01>-|10>)/(√2))((|01>-|10>)/(√2)) =((|0101>-|0110>-|1001>+|1010>)/2) *)
Ψ=KroneckerProduct[EPR,EPR]
MatrixForm[Ψ]

(({{0}, {0}, {0}, {0}, {0}, {(1/2)}, {-(1/2)}, {0}, {0}, {-(1/2)}, {(1/2)}, {0}, {0}, {0}, {0}, {0}}))

(* Bell Basis *)
Φ1=((1/(√2)))*(({{1}, {0}, {0}, {1}}));  (* Φ1 = ((|00>+|11>)/(√2)) *) MatrixForm[Φ1]
Φ2=((1/(√2)))*(({{0}, {1}, {1}, {0}})) ; (* Φ1 = ((|001>+|10>)/(√2)) *) MatrixForm[Φ2]
Φ3=((1/(√2)))*(({{0}, {1}, {-1}, {0}}));  (* Φ1 = ((|001>-|10>)/(√2)) *) MatrixForm[Φ3]
Φ4=((1/(√2)))*(({{1}, {0}, {0}, {-1}}));  (* Φ1 = ((|00>-|11>)/(√2)) *) MatrixForm[Φ1]

(* Bell Projectors *)
P1=Φ1.Φ1\[ConjugateTranspose]; MatrixForm[P1]
P2=Φ2.Φ2\[ConjugateTranspose]; MatrixForm[P2]
P3=Φ3.Φ3\[ConjugateTranspose]; MatrixForm[P3]
P4=Φ4.Φ4\[ConjugateTranspose]; MatrixForm[P4]

(* Testing that these are actually projectors *)
P1.P1==P1
P2.P2==P2
P3.P3==P3
P4.P4==P4

(* Testing for orthogonality *)

P1.P2

(* Bell Measurement of middle 2 qubits of Ψ *)

I2=(({{1, 0}, {0, 1}}))
PP1=KroneckerProduct[KroneckerProduct[I2,P1],I2];MatrixForm[PP1]
PP2=KroneckerProduct[KroneckerProduct[I2,P2],I2];MatrixForm[PP2]
PP3=KroneckerProduct[KroneckerProduct[I2,P3],I2];MatrixForm[PP3]
PP4=KroneckerProduct[KroneckerProduct[I2,P4],I2];MatrixForm[PP4]

{{Prob1}}=Ψ\[ConjugateTranspose].PP1.Ψ; Prob1
Ψ1=((PP1.Ψ)/(√Prob1)); MatrixForm[Ψ1]

(({{-(1/2)}, {0}, {0}, {0}, {0}, {0}, {-(1/2)}, {0}, {0}, {-(1/2)}, {0}, {0}, {0}, {0}, {0}, {-(1/2)}})) (* = -((|0000>+|0110>+|1001>+|1111>)/2) *)

{{Prob2}}=Ψ\[ConjugateTranspose].PP2.Ψ; Prob2
Ψ2=((PP2.Ψ)/(√Prob2)); MatrixForm[Ψ2]

(({{0}, {0}, {0}, {(1/2)}, {0}, {(1/2)}, {0}, {0}, {0}, {0}, {(1/2)}, {0}, {(1/2)}, {0}, {0}, {0}})) (* = ((|0011>+|0101>+|1010>+|1100>)/2) *)

{{Prob3}}=Ψ\[ConjugateTranspose].PP3.Ψ; Prob3
Ψ3=((PP3.Ψ)/(√Prob3)); MatrixForm[Ψ3]

(({{0}, {0}, {0}, {-(1/2)}, {0}, {(1/2)}, {0}, {0}, {0}, {0}, {(1/2)}, {0}, {-(1/2)}, {0}, {0}, {0}})) (* = ((-|0011>+|0101>+|1010>-|1100>)/2) *)

{{Prob4}}=Ψ\[ConjugateTranspose].PP4.Ψ; Prob4
Ψ4=((PP4.Ψ)/(√Prob4)); MatrixForm[Ψ4]

(({{(1/2)}, {0}, {0}, {0}, {0}, {0}, {-(1/2)}, {0}, {0}, {-(1/2)}, {0}, {0}, {0}, {0}, {0}, {(1/2)}})) (* = ((|0000>-|0110>-|1001>+|1111>)/2) *)
Ψ3=((PP3.Ψ)/(√Prob3)); MatrixForm[Ψ3]

({{0, Resulting State, Probability}, {Measure wrt P1, Result1=-((|0000>+|0110>+|1001>+|1111>)/2), (1/4)}, {Measure wrt P2, Result2=((|0011>+|0101>+|1010>+|1100>)/2), (1/4)}, {Measure wrt P3, Result3=((-|0011>+|0101>+|1010>-|1100>)/2), (1/4)}, {Measure wrt P4, Result4=((|0000>-|0110>-|1001>+|1111>)/2), (1/4)}})

The renaming of the qubits given by the permutation (({{0, 1, 2, 3}, {1, 2, 3, 0}}))=(0123000000) induces the permutation of the basis states given by 
                 |b3,b2,b1,b0> --> |b0,b3,b2,b1>
This induced permutation is  Perm =(0)(1,2,48)(5,10)(6,12,9,3)(7,14,13,11)(15), 
where, for example, 9 stands for 1001.After applying Prm, we have:

({{0, Resulting State, Probability}, {Measure wrt P1, Result1'=-((|0000>+|1100>+|0011>+|1111>)/2)=-((|00>+|11>)/(√2)).((|00>+|11>)/(√2)), (1/4)}, {Measure wrt P2, Result2'=((|0110>+|1010>+|0101>+|1001>)/2)=((|01>+|10>)/(√2)).((|01>+|10>)/(√2)), (1/4)}, {Measure wrt P3, Result3'=((-|0110>+|1010>+|0101>-|1001>)/2)=((|01>-|10>)/(√2)).((|01>-|10>)/(√2)), (1/4)}, {Measure wrt P4, Result4'=((|0000>-|1100>-|0011>+|1111>)/2)=((|00>-|11>)/(√2)).((|00>-|11>)/(√2)), (1/4)}})

Thus,before the, qubits 0 & 1  and qubits 2 & 3 respectivly form EPR pairs.
But after the measurement, qubits 0 & 2 and qubits 1 & 3 respectively form EPR pairs.
Hence we have "swapped" the entanglement.