## 3 Construction of Algorithm  
In this section, we construct an algorithm to get the approximate solution of the Bloch equation. Using this algorithm, we can then obtain magnetisation in each direction. Let us take the following approximations:  
\[\frac{d^{x}N_{x}(t)}{dt^{\alpha}}=C_{1}^{T}q(t),\ \frac{d^{\beta}N_{y}(t)}{dt^{ \beta}}=C_{2}^{T}q(t),\ \frac{d^{\gamma}N_{z}(t)}{dt^{\gamma}}=C_{3}^{T}q(t), \tag{10}\]  
and  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]From Equations (10) and (11), we can write  
\[N_{x}(t)=C_{1}^{T}I^{(\alpha)}q(t)+L^{T}q(t), \tag{12}\]  
\[N_{y}(t)=C_{2}^{T}I^{(\beta)}q(t)+M^{T}q(t), \tag{13}\]  
\[N_{z}(t)=C_{3}^{T}I^{(\gamma)}q(t)+N^{T}q(t) \tag{14}\]  
Using Equations (10), (12), (13) and (14) in Equation (3), we get  
\[C_{1}^{T}\Big{(}I+\frac{1}{T_{2}}I^{(\alpha)}\Big{)}-w_{0}C_{2}^{T}I^{(\beta)} =w_{0}M^{T}-\frac{1}{T_{2}}L^{T} \tag{15}\]  
\[w_{0}C_{1}^{T}I^{(\alpha)}+C_{2}^{T}\Big{(}I+\frac{1}{T_{2}}I^{(\beta)}\Big{)} =-w_{0}L^{T}-\frac{1}{T_{2}}M^{T} \tag{16}\]  
\[C_{3}^{T}\Big{(}I+\frac{1}{T_{1}}I^{(\gamma)}\Big{)}=O^{T}-\frac{1}{T_{1}}N^{T} \tag{17}\]  
where \(I^{(\alpha)}\), \(I^{(\beta)}\) and \(I^{(\gamma)}\) are the operational matrices of non-integer-order integration of order \(\alpha,\beta\) and \(\gamma\), respectively. Here \(I\) is an identity matrix. The simpler form for Equations (15)-(17) is given as follows:  
\[C_{1}^{T}W_{1}-C_{2}^{T}W_{5}=G_{1}, \tag{18}\]  
\[C_{1}^{T}W_{4}+C_{2}^{T}W_{2}=G_{2}, \tag{19}\]  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]0  
where  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]1  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]2  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]3  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]4  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]5  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]6  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]7  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]8  
The matrices \(W_{1},\ W_{2},\ W_{3},W_{4},\ W_{5},\ G_{1},\ G_{2}\) and \(G_{3}\) are given in terms of known values, so these matrices are known matrices. On solving Equations (18)-(20), we get  
\[N_{x}(0)=L^{T}q(t),\ N_{y}(0)=M^{T}q(t),\ N_{z}(0)=N^{T}q(t),\ \frac{N_{0}}{T_{1}}=O^{T}q(t). \tag{11}\]9  
\[N_{x}(t)=C_{1}^{T}I^{(\alpha)}q(t)+L^{T}q(t), \tag{12}\]0  
\[N_{x}(t)=C_{1}^{T}I^{(\alpha)}q(t)+L^{T}q(t), \tag{12}\]1  
In Figures 3 and 4, we have presented 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order \((\alpha=\beta=\gamma=0.9)\), respectively. These figures show the dynamics of \(N_{x},\ N_{y}\) and \(N_{z}\) for fractional-order relaxation. In Figure 3, the entire trajectory of magnetisation is shown in 3D for fractional order \((\alpha=\beta=\gamma=0.9)\) starting at i. c. \(\left(N_{x}(0),N_{y}(0),N_{z}(0)\right)\) and returning to \(N_{0}\). From Figure 4, it is clear that the magnetisation \(N_{x}\) in \(x-\)direction increases with time, and the magnetisation \(N_{y}\) in \(y-\)direction decreases with time. In Figures 5 and 6, we have presented the 3D and 2D plots of the numerical solutions of the Bloch equation for fractional order (\(\alpha=\beta=\gamma=0.8\)), respectively. These figures show the dynamics of \(N_{x},\ \ N_{y}\) and \(N_{z}\) for the fractional-order relaxation.