spin = 3/2; MagnusOrder = 0;
(* parameters of the non-perturbed interaction *)
wL = 2 * Pi * 79.46 * 10^6; thetaZ = 0; phiZ = 0;
alphaQ = 0; betaQ = 0; gammaQ = 0; eta = 0; cQ = 2 * (2 * Pi * 51.0 * 10^3);
H0 = HQ[spin] + HZ[spin];

(* parameters of the RF interaction *)
wRF = wL; w1 = Pi/(2*10^-6); thetaRF = Pi/2; phiRF = 0; phase = Pi/2;
H1 = HRF[spin];

(* matrices of fictitious spin 1/2 operators *)
I12 = {{0, 1/2, 0, 0}, {1/2, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
I23 = {{0, 0, 0, 0}, {0, 0, 1/2, 0}, {0, 1/2, 0, 0}, {0, 0, 0, 0}};
I34 = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1/2}, {0, 0, 1/2, 0}};

(* determining the initial state *)
rho0 = {{1, 0, 0, 0}, {0, 3, 0, 0}, {0, 0, -3, 0}, {0, 0, 0, -1}};

(* calculation of relative intensity versus flip angle *)
NumPoints = 30; PlotMax = 90;
M12 = Table[{0, 0}, {NumPoints + 1}];
M23 = Table[{0, 0}, {NumPoints + 1}]; M34 = Table[{0, 0}, {NumPoints + 1}];
Do[flip = i * (PlotMax/NumPoints); TP = (flip/90) * (2 * 10^-6);
rho1 = Pulse[rho0];
j12 = ((MatrixExp[(I*H0*TP)/hbar].I12.MatrixExp[(- I*H0*TP)/hbar]));
j23 = ((MatrixExp[(I*H0*TP)/hbar].I23.MatrixExp[(- I*H0*TP)/hbar]));
j34 = ((MatrixExp[(I*H0*TP)/hbar].I34.MatrixExp[(- I*H0*TP)/hbar]));
M12[[i + 1]] = {flip, Sqrt[3]*Tr[rho1.j12]};
M23[[i + 1]] = {flip, 2*Tr[rho1.j23]};
M34[[i + 1]] = {flip, Sqrt[3]*Tr[rho1.j34]},
{i, 1, NumPoints}]

(* plot results *)
ListPlot[{M12, M23, M34}, PlotRange -> {-1, 5},
TicksStyle -> Directive[20], PlotJoined -> True, PlotStyle -> Thick]