Accurate Bloch and Bloch-McConnell model simulation is essential for proper modelling and becomes crucial for iterative RF pulse design methods. This work demonstrates the advantages of a piecewise analytical solution using eigenvalues and eigenvectors and an asymmetric and symmetric operator splitting scheme. Those three solvers are compared to five other numerical solution methods regarding accuracy and runtime. The symmetric operator splitting prevails in both compared to the others. Moreover, it is quadratic convergent with respect to the time step and therefore seems to be a good choice for optimal control approaches for Bloch and for Bloch-McConnell models.