MATLAB: Simple For cycle (a small problem with the indexes)

forfor loopMATLAB

I have a small problem with the last loop in the code:
clear
clc
close all
load matlab_10.mat
init= 24500;
%% Gas properties
k = 1.4;          %Gas constant ratio, cp/cv()
R = 287.05;       %Ideal gas constant [J/Kg*K]
p01 = 1e5;        %Compressor inlet pressure(Pa)
T01 = 303.35;     %Compressor inlet temperature(K)
cp = 1005;        %Specific heat capacity for constant pressure, property of gas(J/(kg*K))
ro1 = 1.15;       %Gas density(kg/m^3)
AFR = 14.71;      %Air Fuel Ratio
%% Turbine's parameters
A_eff = .03;       %[m^2]
A_0 = A_eff/.6;
E_R = [1:0.001:4]; %Expansion ratio
T_3 = 1173.15;     %Discharge temperature [K]
eta_t = 0.6;
%% Compressor's parameters
eta_c = 0.95;
I = 0.001;                                  %Moment of inertia for the compressor spool(kg*m^2)
mass_flow_compressor = 0.273675009794891;   %adimensional mass flow from Greitzer steady-state [kg/s]
% m_c = mass_flow_compressor*(ro1*Ac*U);    %actual compressor's mass flow [kg/s]  
mass_flow_compressor_corrected = (mass_flow_compressor*(T01^0.5))./p01; %[m*s*(k)^0.5]

P_R = 2.591609007038750; %pressure ratio from Greeitzer
p2 = P_R*p01;  %[Pa] 
%% 
mass_flow_compressor_corrected = (mass_flow_compressor_corrected*ones(1,length(E_R))).*(1+(1./AFR));
mass_flow_turbine_corrected = A_eff.*((k./R).^0.5).*((1./E_R).^(1/k)).*(((2./k-1).*(1-(1./E_R).^((k-1)./k))).*0.5);
% mass_flow_corrected = (mass_flow).*p01;
x_ER = interp1(mass_flow_turbine_corrected - mass_flow_compressor_corrected, E_R, 0);
y_ER = A_eff.*((k./R).^0.5).*((1./x_ER).^(1/k)).*(((2./k-1).*(1-(1./x_ER).^((k-1)./k))).*0.5);
figure()
plot(E_R,mass_flow_turbine_corrected,...
    E_R,mass_flow_compressor_corrected,'linewidth',3)
grid on
grid minor
xlabel('ER')
legend('turbine','compressor')
ylabel('m*s*K^{0.5}')
title('ER vs mass flow')
hold on
plot(x_ER, y_ER, 'pg','linewidth',3)
hold off
text(x_ER, y_ER, sprintf('\\uparrow\nIntersection:\nx\\_ER = %.3f\ny\\_ER = %.3f', x_ER, y_ER),...
    'HorizontalAlignment','left', 'VerticalAlignment','top')
%% Power Steady Case
Pt = ((y_ER.*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER.^((k-1)./k))));
Pc = ((y_ER.*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
%% Pulsating
mass_flow_compressor_vet = phi_10(init:end);
mass_flow_compressor_vet_corrected = (mass_flow_compressor_vet*(T01^0.5))./p01; %[m*s*(k)^0.5]
mass_flow_compressor_vet_corrected = mass_flow_compressor_vet_corrected*ones(1,length(E_R)).*(1+(1./AFR));
for i=1:length(phi_10(init:end))
%        
   x_ER_vet(i,:) = interp1(mass_flow_turbine_corrected - mass_flow_compressor_vet_corrected(i,:), E_R, 0);
   y_ER_vet(i,:) = A_eff.*((k./R).^0.5).*((1./x_ER_vet(i,:)).^(1/k)).*(((2./k-1).*(1-(1./x_ER_vet(i,:)).^((k-1)./k))).*0.5);
% 
%    Pt_vet = ((y_ER_vet.*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER_vet.^((k-1)./k))));
%    Pc_vet = ((y_ER_vet.*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
%    omega_vet = (2.*((Pt_vet-Pc_vet)./(I))+55000.^2).^0.5;
end
figure()
plot(t_10(init:end),omega_vet,'linewidth',3)
grid on
grid minor
xlabel('t [s]')
legend('\omega')
ylabel('rpm')
hold on
title('\omega')
my problem is I have the first omega, let's call it omega(1), and by that iteration I need to find the others.
The algorithm I wrote on paper foor the first two steps is in the following figure: (anyway I've already written the code in the comment inside the for-loop, i just miss the indexes)

Best Answer

There you go: 
omega_vet(1) = 55000;
for i=1:length(phi_10(init:end))
       
   x_ER_vet(i) = interp1(mass_flow_turbine_corrected - mass_flow_compressor_vet_corrected(i), E_R, 0);
   y_ER_vet(i) = A_eff.*((k./R).^0.5).*((1./x_ER_vet(i)).^(1/k)).*(((2./k-1).*(1-(1./x_ER_vet(i)).^((k-1)./k))).*0.5);
   Pt_vet(i)   = ((y_ER_vet(i).*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER_vet(i).^((k-1)./k))));
   Pc_vet(i)   = ((y_ER_vet(i).*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
   
   omega_vet(i+1) = (2.*((Pt_vet(i)-Pc_vet(i))./(I)) + omega_vet(i).^2).^0.5;
end
figure()
plot(t_10(init:end),omega_vet(1:end-1),'linewidth',3)
Related Question