Each step in ReferenceModel has to use an intiial state that is the final state from the previous step. Using that term y0 is incorrect. In order to keep track of the state vector, it appears a ss object is required as the input in the call to lsim()
s = tf('s');
Gm = 2163/(s^2 + 80*s + 2163);
samplingTime = 2e-3;
totalTime = 0.5;
time = 0:samplingTime:totalTime;
r = 80;
% Using for loop
ym(1,:) = 0; xm = [0;0]; % initial state
for k = 2:length(time)
[ym(k,:),xm] = ReferenceModel(Gm,r,samplingTime,xm); % send in state from prvevious time step
end
% Without for loop
ym2 = lsim(Gm,r*ones(size(time)),time);
hold on
plot(time,ym,'bo')
plot(time,ym2,'r--')
legend('ym','ym2')
function [ym,xm] = ReferenceModel(Gm,r,samplingTime,x0)
tspan = 0:samplingTime/10:samplingTime;
[ym,~,xm] = lsim(ss(Gm),r*ones(size(tspan)),tspan,x0); % use ss form and use initial state correctly
ym = ym(end);
xm = xm(end,:);
end
