以红色标示一下你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, '