In this paper we are concerned with the problem of solving numerically isospectral flows. These flows are characterized by the differential equation
\[
L
′
=
[
B
(
L
)
,
L
]
,
L
(
0
)
=
L
0
,
L’ = [B(L), L], \quad L(0)=L_0,
\]
where
L
0
L_0
is a
d
×
d
d\times d
symmetric matrix,
B
(
L
)
B(L)
is a skew-symmetric matrix function of
L
L
and
[
B
,
L
]
[B,L]
is the Lie bracket operator. We show that standard Runge–Kutta schemes fail in recovering the main qualitative feature of these flows, that is isospectrality, since they cannot recover arbitrary cubic conservation laws. This failure motivates us to introduce an alternative approach and establish a framework for generation of isospectral methods of arbitrarily high order.