各位大神，求PID算法的matlab程序，救急啊！

clear
%
clear u y
n=5;
Nmax=300;

for k=1:Nmax,
u(k)=1; y(k)=0;
end
%
k=0;
e=0;
sum=0;
% loop:
for i=n+1:Nmax,
% The (actually unknown!) process, for example. Feel free to change this
% as you please.
Ts=1; % Sample time (seconds)
Kp=3; % The process gain
T=20-10*(i/2000); % The first time-constant(seconds)
ksi=0.05+0.05*(i/2000);
wn=1;%

% Translating the system to the discrete-time
num=Kp;
den=conv([T 1],[1/wn^2 2*ksi/wn 1]);
[B,A]=c2dm(num,den,Ts);
%
% u(i)=5*sign(u(i-1));
% k=k+1;
% if k>50,
% k=0;
% u(i)=-5*sign(u(i-1));
%end

% The process
y(i)=-A(2:length(A))*y(i-1:-1:i-length(A)+1)'+B(2:length(B))*u(i-1:-1:i-length(B)+1)';
% the controller
Kc=0.26;
Ti=6;
Td=0.1;
setpoint=1;
eold=e;
e=setpoint-y(i);
sum=sum+e;
u(i)=Kc*(e+Ts/Ti*sum+Td/Ts*(e-eold));
%
end
figure(1)
subplot(211), plot(u),ylabel('input'),subplot(212),plot(y),ylabel('output')