Hi @Sneh
That's because the response of one of the states, 
exploded. This it cannot run to the end 
sec.
exploded. This it cannot run to the end sec.% t = 0:0.01:(2*pi*4);
t = 0:0.1:3.9;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[t, q] = ode45(@euler_eqns, t, [0.5; 0; 0; 0.866; 0; -1; 0]);
% [t, q] = ode45(@euler_eqns, t, [0.5; 0; 0; 0.866; 0; -1; 0], options);
plot(t, q(:,5))
% K = q(:,1).^2 + q(:,2).^2 + q(:,3).^2 + q(:,4).^2;
% K0 = 1;
% c = K - K0;
% plot(t, K);
% xlabel('Time (s)');
% ylabel('K');
% plot(t, c);
% xlabel('Time (s)');
% ylabel('K-K0');
function dqwdt = euler_eqns(t,qwd, q)
q = qwd(1:4);
w = qwd(5:7);
I = 500;
J = 125;
K = (I-J)/I;
T = 35;
ohm = 1;
qstar = 1/ohm;
dqdt = zeros(4,1);
dqdt(1) = pi*( q(4)*w(1) + q(3)*(qstar - w(2) - ohm) + q(2)*w(3));
dqdt(2) = pi*( q(3)*w(1) + q(4)*(w(2) - ohm - qstar) - q(1)*w(3));
dqdt(3) = pi*(-q(2)*w(1) - q(1)*(qstar - w(2) - ohm) + q(4)*w(3));
dqdt(4) = pi*(-q(1)*w(1) - q(2)*(w(2) - ohm - qstar) - q(3)*w(3));
dwdt = zeros(3,1);
dwdt(1) = 12*pi*K*(q(2)*q(3) + q(1)*q(4))*(1 - 2*q(1)^2 - 2*q(2)^2) - 2*pi*K*w(2)*w(3) + 2*pi*qstar*w(3);
dwdt(2) = 0;
dwdt(3) = -24*pi*K*(q(3)*q(1) - q(2)*q(4))*(q(2)*q(3) + q(1)*q(4)) + 2*pi*K*w(1)*w(2) + 2*pi*qstar*w(1);
% Combine the time derivatives into a single vector
dqwdt = [dqdt;
dwdt];
end
