hi ,i have this function that ode45 integrate, the numerical integration is correct, but in the main script I have to plot the vector 'alfa', and in order to do this i have to rewrite the same code, but i need to do a cicle 'for', how is this cicle ?

function [ derivate ] = eq_moto_vincolo_q_piatto(t, X )

format long
global mi m CQ T w A rho0 H R_earth beta a0 a1 b0 b1 b2 c0 c1 c2 c3  
global tD_ott alfaD_ott sigmaD_ott
global t_pre alfa_pre

%vettore incognita

r=X(1);
lambda=X(2);
lat=X(3);
gamma=X(4);
v=X(5);
zita=X(6);
calore_assorbito=X(7);

%%equazioni del moto

%equazioni cinematiche

if r<R_earth
  v=0;
end

dr=v*sin(gamma);
dlambda=v*cos(gamma)*cos(zita)/(r*cos(lat));
dlat=v*cos(gamma)*sin(zita)/r;

%densità

rho=rho0*exp(-(r-R_earth)/H);

%parametri aerodinamici

alfa_tilde=interp1(tD_ott,alfaD_ott,t,'spline');

flusso_tilde = ( c0 + c1.*(alfa_tilde) + c2.*(alfa_tilde).^2 + c3.*(alfa_tilde).^3 ) * ( 17700 .* sqrt(rho * 0.002970/1.75e9 ) .* (0.0001.*v * 3280.84 ).^(3.07) );    %BTU/ft^2 sec  

if flusso_tilde < 70
    alfa = alfa_tilde
else

      %c0 + c1*(alfa) + c2*(alfa)^2 + c3*(alfa)^3 = termine_noto ;

      termine_noto = 70 / ( 17700 * sqrt(rho * 0.002970/1.75e9 ) * (0.0001 * v * 3280.84 )^(3.07) ) ;
      coefficienti = [c3 c2 c1 c0-termine_noto];

alfa_sol = roots(coefficienti)

ind_alfa_sol=find(imag(alfa_sol));
alfa_sol_complessi=alfa_sol(ind_alfa_sol);
alfa_sol(ind_alfa_sol)=[]

ind = find ( min(abs(alfa_sol - alfa_pre)) );
alfa = alfa_sol(ind)

end

sigma_int=interp1(tD_ott,sigmaD_ott,t,'spline');
sigma = sigma_int * pi/180;

CL = a0 + a1*(alfa);
CD = b0 + b1*(alfa) + b2*(alfa)^2;

%forza aerodinamica

L=0.5*rho*CL*A*v^2;
D=0.5*rho*CD*A*v^2;
Q=0.5*rho*CQ*A*v^2;

%equazioni dinamiche

dgamma= ( -mi/(r^2*v) + v/r ) * cos(gamma) + ( T*cos(beta)*sin(alfa)+L*cos(sigma)+Q*sin(sigma) )/(m*v) + (w^2*r/v)* cos(lat)*( cos(lat)*cos(gamma) + sin(lat)*sin(gamma)*sin(zita) ) +2*w*cos(lat)*cos(zita); 
dv= -(mi/r^2)*sin(gamma) + ( T*cos(beta)*cos(alfa)-D )/(m) + w^2*r*cos(lat) * ( cos(lat)*sin(gamma) - sin(lat)*cos(gamma)*sin(zita) )  ; 
dzita= -v/r * tan(lat)*cos(gamma)*cos(zita) + ( T*sin(beta)-L*sin(sigma)+Q*cos(sigma) )/( m*v*cos(gamma) ) + 2*w*cos(lat)*tan(gamma)*sin(zita) - (w^2*r)/(v*cos(gamma))*sin(lat)*cos(lat)*cos(zita) - 2*w*sin(lat); 

%flusso di calore

flusso = ( c0 + c1.*(alfa) + c2.*(alfa).^2 + c3.*(alfa).^3 ) * ( 17700 .* sqrt(rho * 0.002970/1.75e9 ) .* (0.0001.*v * 3280.84 ).^(3.07) )   %BTU/ft^2 sec  

if t>t_pre
    t_pre = t;
    alfa_pre = alfa;
end

%funzione da integrare

derivate=[dr; dlambda; dlat; dgamma; dv; dzita; flusso];

[Walter Roberson]

You can use a global or a (better) a shared variable to keep a record of alfa, probably updating at your "if t>t_pre" . You should probably keep a record of the associated t value as well.