EDIT: the OP has provided the answer in a comment, see the bottom of this post. Most of my answer (until the dividing line) is irrelevant to
resolving the problem, but I'm leaving it in place as it provides
background information and advice on checking the numerical method.
I'm going to try and interpret your equations and offer an explanation. Your notation is different from what I am used to, so I can't guarantee that we'll agree on everything. I'm working from H Goldstein "Classical Mechanics".
Your first two equations together express the relation between angular momentum $\vec{K}$ and applied torque $\vec{M}$ in the stationary, non-rotating, space-fixed frame. So, you seem to be using the subscript $G$ to denote this frame,
or perhaps it is to indicate moments taken about the centre of mass.
To avoid ambiguity, I am going to use $S$ for this space-fixed frame.
$$
\frac{d \vec{K}_S}{d t} = \vec{M}_S
$$
This is the basis of your impulse-momentum formula, which applies in the stationary frame. It is possible to integrate this equation forward
$$
\vec{K}_S(t+\Delta t) \approx \vec{K}_S(t) + \vec{M}_S\Delta t
$$
Some people do tackle the rotational equations of motion this way.
But as I think you have noticed, obtaining $\vec{\omega}_S(t+\Delta t)$
involves taking account of the rotation of the inertia tensor (in space-fixed coordinates), over the timestep. You haven't indicated how you are
handling this part of the equations of motion.
Certainly one can write
$$
\vec{\omega}_S(t)=I_S^{-1}(t)\vec{K}_S(t),
\qquad
\vec{\omega}_S(t+\Delta t)
=I_S^{-1}(t+\Delta t)
\vec{K}_S(t+\Delta t)
$$
if one has $I_S$ at both times.
When you write
$$
\vec{K}_{G_1} - \vec{K}_{G_0} = I(\vec{\omega}_1-\vec{\omega}_0)
$$
it is not clear how you are transforming the body-fixed quantities on the right,
with $I$ assumed constant,
into the space-fixed form needed on the left,
or indeed what $G_0$ and $G_1$ are.
So I have some doubts about these equations.
The time derivative of the angular momentum in rotating, body-fixed coordinates
may be written (following Goldstein, in your notation, but using the subscript $B$ to indicate vectors resolved in the body-fixed frame)
$$
\frac{d \vec{K}_B}{d t} + \vec{\omega_B}\times\vec{K}_B = \vec{M}_B
$$
For these quantities
we can write $\vec{K}_B=I_B\vec{\omega}_B$ with $I_B$ independent of time.
This is your second equation.
It is possible to solve the equations of motion in this form,
rearranging them as you have done
$$
\frac{d \vec{\omega}_B}{d t}
= I_B^{-1}\left[\vec{M}_B -\vec{\omega_B}\times I_B\vec{\omega}_B\right]
$$
I believe that these equations are correct.
That leads me on to the question:
how would you test whether it is working properly?
I'm not sure how important this is for your application.
In quantitative physics applications,
to test the solutions, it is helpful to check some things that are meant to be conserved, rather than just relying on the visual appearance of the dynamics.
Suppose we consider one rotating body, being acted on by a torque $\vec{M}$
which is derived from a potential energy function. Then the equations of motion should conserve the total energy (kinetic + potential, or KE+PE). This can be seen because
$$
\frac{d\, \text{KE}}{dt} = \vec{\omega}_B\cdot\vec{M}_B =\vec{\omega}_S\cdot\vec{M}_S
$$
and it can be shown that the right hand side is the negative of the rate of
change of PE, because of the way $\vec{M}$, the torque, is defined.
However, if the body-fixed equations of motion are solved
without the $\vec{\omega_B}\times I_B\vec{\omega}_B$ term,
they will still conserve energy,
because this term is perpendicular to the angular velocity.
It may be worth checking how well your method conserves
kinetic energy in the special case of zero torque.
However,
I believe that the incorrect equations (missing that term)
will not satisfy conservation of total angular momentum $\vec{K}_S$
in the space-fixed frame.
It may be possible to test this with zero external torque $\vec{M}_S$
on the body.
For a general inertia tensor $I$,
so $\vec{\omega}$ and $\vec{K}$ are not parallel,
and $\vec{M}_S=0$
you should see the body wobble around, so $\vec{\omega}$ changes noticeably
in the space-fixed frame,
but $\vec{K}_S$ should be constant.
I recommend doing that test, as well as the energy conservation test.
Just to recap: I'm not sure that your second method
based on the impulse-momentum formula,
is just the body-fixed equation missing the $\vec{\omega_B}\times I_B\vec{\omega}_B$ term; it may be something more complicated than that.
If so, my comments about the consequences of omitting that term may be irrelevant.
It is certainly possible to solve the rotational equations,
by applying the impulse-momentum form in space-fixed coordinates,
and taking care with the rotation of the inertia tensor.
But I do have my doubts about the way this is implemented.
And it is also possible to solve the equations in body-fixed form,
as you seem to be doing,
correctly as far as I can judge.
In this case, my only worries are about whether the numerical method
is accurate enough, conserving energy and momentum.
For this form of the equations (with $\vec{\omega}$ appearing
on the right) it may be that predictor-corrector methods are
more accurate than explicit Euler.
And then, I have no idea how you are tackling the
associated orientational equations of motion,
which actually rotate your rigid body.
Of course these details may not be critical in your application.
EDIT following OP comment.
In the event that the applied torque takes the form of an
instantaneous impulse $\vec{M}=\vec{Q}\delta(t)$, where $\vec{Q}$ is a constant vector, and $\delta(t)$ is a Dirac delta function,
then the apparent discrepancy disappears.
The OP has correctly resolved the problem in this case.
In the space-fixed frame, time integrating over the infinitesimal interval
during which the impulse is applied gives
$$
\vec{K}_S' - \vec{K}_S = \vec{Q}_S \qquad
\Rightarrow \qquad
\vec{\omega}_S' - \vec{\omega}_S = I_S^{-1}\vec{Q}_S
$$
where $'$ denotes the values after the impulse.
The space-fixed inertia tensor does not change during this infinitesimal
time interval.
In the body-fixed frame the equation for $d\vec{K}_B/dt$ includes an extra term $\vec{\omega_B}\times\vec{K}_B$. However, during the infinitesimal interval of the impulse, $\vec{\omega_B}$ and $\vec{K}_B$ will have only
step function discontinuities, no delta-functions, and so integrating them over this interval will give zero. Therefore
$$
\vec{K}_B' - \vec{K}_B = \vec{Q}_B \qquad
\Rightarrow \qquad
\vec{\omega}_B' - \vec{\omega}_B = I_B^{-1}\vec{Q}_B
$$
where everything is expressed in body-fixed axes.
So, for a delta-function impulse,
no error arises from using this equation.
It is easy to check angular momentum conservation (in the space-fixed frame) by converting these body-fixed equations to that frame,
using relations such as $\vec{Q}_B= R\,\vec{Q}_S$ and $I_B=R\, I_S\, R^T$,
where $R$ is the rotation matrix describing the transformation.
$or$$, like so $ \frac{dKg}{dt} = I * \frac{dw}{dt} + w * Iw $ $\endgroup$