for k=1:4
S(k)=(1/(k+1))*(m-b*y-m*r); % How to get y here?
I(k)=(1/(k+1))*(a*e-(m+r)*i);
end
S=sort(S);
I=sort(I,'descend');
%Now Calculate Y
Or
y=zeros(1,4)
for k=1:4
S(k)=(1/(k+1))*(m-b*y(k)-m*r);
I(k)=(1/(k+1))*(a*e-(m+r)*i);
y(k+1)=y(k)+(S(k)*I(k));
end
S=sort(S);
I=sort(I,'descend');
% Canculate y again % I didnot find any sense here
