以紅色標示一下你code有問題的地方(我修正的地方)
clear
Ts=25;
p=1.00151*10^(-4);
q=8.67973*10^(-3);
r=40.0;
a=exp(-p*Ts);
b=(q/p)*(1-a);
yo =25;
% 以下幾行只是幫你修改成易讀的code
ref = [35*ones(1, 40), 50*ones(1, 40), 65*ones(1, 40), 80*ones(1, 380)];
GE=1/15;
GC=1/15;
GU=450;
y = zeros(1, 500);
y(1)=25;
e = zeros(1, 500);
ee = zeros(1, 500);
ce = zeros(1, 500);
ce(1)=0;
cee = zeros(1, 500);
mu = zeros(1, 500);
u = zeros(1, 500);
%
for k=1:500
e(k)=ref(k)-y(k);
ee(k)=e(k)*GE;
if k >= 2
ce(k)=e(k)-e(k-1);
end
cee(k)=e(k)*GC;
mu(k)=1;
u(k)=mu(k)*GU;
y(k+1)=a*y(k)+b/(1+exp(0.5*y(k)-r))*u(k)+(1-a)*yo;
end
% 你的y長度是501，其他參數e, ee, ce, cee, mu, u的長度都是500
% 所以0:25:499 * Ts長度是500才能畫，y只取前五百個畫
subplot(2,2,1)
plot(0:25:499 * Ts, y(1:500), 0:25:499 * Ts, ref, '