这是一个别人写的z—n法整定pid的程序，求大神指点指点

%%%%Program
function [Kp,Ti,Td]=Z_N_PIDparameter(num,den,delay,t_start,t_end)
h=.0001;
t=t_start:h:t_end;                     
if delay==0                             %no delay               
        sys=tf(num,den);
else
    sys=tf(num,den,'iodelay',delay);    %containing delay        
end                                    
k=dcgain(sys);                          %get the gain of the system                                                                 
y=step(sys,t);
dy=diff(y)/h;                                  %derivation         
[m,n]=max(dy);                          
t_col=size(t,2);                        
t_max=(t_end-t_start)*n/t_col;
y_index=t_col*t_max/(t_end-t_start);
y_tangentline=y(round(y_index));
b_intercept=y_tangentline-m*t_max;        %b_intercept
syms t;
tk=solve(k-m*t-b_intercept);
tk=double(tk);
L=t_max;
T=tk-t_max;
alafa=k*L/T;
Kp=1.2/alafa;
Ti=2*L;
Td=L/2;
%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%Function Help
%%%%%%%%%%Example 1
%
%Transfer function(containing timedelay):
%                        20
%   exp(-0.5*s) * -------------------
%                    1.6 s^2 + 4.4 s + 1
%
%
%[Kp,Ti,Td]=Z_N_PIDparameter(20,[1.6 4.4 1],0.5,0,30)
%we can get
%                      Kp = 0.1733
%               Ti =3.0468
%               Td =0.7617
%


%%%%%%example 2
%Transfer function(no delay):
%
%           20
%   -------------------
%   1.6 s^2 + 4.4 s + 1
%
%[Kp,Ti,Td]=Z_N_PIDparameter(20,[1.6 4.4 1],0,0,30)
%we can get
%                Kp = 0.2580
%                Ti =2.0468
%                Td =0.5117
%
对y求导的时候dy=diff(y)/h,后面为什么要除以h？
y_index=t_col*t_max/(t_end-t_start);还有这一步，这个y_index是切点的y坐标吗？
后面变量好像也有点问题。
来点支援！！！！