Extracting an extra Matrix from ode45 solver

Hey guys, I want to extract an extra matrix (eta_d) from the solver, which will contain values calculated at each time step.. Is there anyway I can do that? Thanks in advance for help.

Here is the function and executive file:

function dydt = ODEsystem(t,y)
phi=sin(t);
eta_d=cos(phi/2);
a=[1/sqrt(3);1/sqrt(3);1/sqrt(3)];
epsilon_d=a*sin(phi/2);
w_d=cos(t)*a;
w_d_skew=[0 -w_d(3,1) w_d(2,1);w_d(3,1) 0 -w_d(1,1);-w_d(2,1) w_d(1,1) 0];
K_p=[30 0 0;0 30 0;0 0 30];
K_d=[45 0 0;0 45 0;0 0 45];
I=[25000 0 0;0 100000 0;0 0 110000];
Iomega_d_dot=I*(-sin(t)*a);           
omegacross_dIsame=w_d_skew*I*w_d;     
dydt = zeros(7,1);
dydt(1) =  (Iomega_d_dot(1,1)+ omegacross_dIsame(1,1)-K_p(1,1)*(eta_d*y(5)+y(6)*epsilon_d(3,1)-epsilon_d(2,1)*y(7)-y(4)*epsilon_d(1,1))-K_d(1,1)*(y(1)-w_d(1,1))+I(2,2)*y(2)*y(3)-I(3,3)*y(2)*y(3))/I(1,1);    
dydt(2) =  (Iomega_d_dot(2,1)+ omegacross_dIsame(2,1)-K_p(2,2)*(eta_d*y(6)-y(5)*epsilon_d(3,1)+epsilon_d(1,1)*y(7)-y(4)*epsilon_d(2,1))-K_d(2,2)*(y(2)-w_d(2,1))-I(1,1)*y(1)*y(3)+I(3,3)*y(3)*y(1))/I(2,2);                                                                               
dydt(3) =  (Iomega_d_dot(3,1)+ omegacross_dIsame(3,1)-K_p(3,3)*(eta_d*y(7)-y(6)*epsilon_d(1,1)+epsilon_d(2,1)*y(5)-y(4)*epsilon_d(3,1))-K_d(3,3)*(y(3)-w_d(3,1))-I(2,2)*y(2)*y(1)+I(1,1)*y(2)*y(1))/I(3,3);          
dydt(4) = -0.5*(y(1)*y(5)+y(2)*y(6)+y(3)*y(7));         
dydt(5) = -0.5*(y(2)*y(7)-y(3)*y(6)-y(4)*y(1));          
dydt(6) = -0.5*(y(3)*y(5)-y(1)*y(7)-y(4)*y(2));
dydt(7) = -0.5*(y(1)*y(6)-y(2)*y(5)-y(4)*y(3));
end

y0 = [0.1;0.3;0.1;1;0;0;0]; 
tspan = [0 20]; 
[t,y] = ode45(@ODEsystem,tspan,y0);

Answer 1

[Torsten]

y0 = [0.1;0.3;0.1;1;0;0;0]; 
tspan = [0 20]; 
[t,y] = ode45(@ODEsystem,tspan,y0);
eta_d = cos(sin(t)/2);

Extracted.

[David Togonidze]

I was trying something similiar but getting the results that were not what i expected. But now I know that it is most likely bug in my code. Thank you for the help.