function K=ljiewding(Gsys,Hsys)
%系統臨界穩定K值的判定，由函數ljiewding和chdihshu組成。Gsys：前向通道傳遞函數；
%Hsys：反饋通道傳遞函數；Gsys、Hsys：TF模型
k=0.0001:1:1000;nn=length(k);        %K取值以1為步長，初步判定臨界不穩定K的值
for i=1:1:800
   p=chdihshu(Gsys,Hsys,k(i));    %當系統增益為k(i)時，求系統特征方程的根向量
   ii=find(real(p)>=0);n=length(ii);
   if(n>0),break,end                                %若找到K值，退出循環
end
if n==0                                        %若沒有找到位于S右半平面的根，則n=0
   disp('系統穩定K的取值范圍為K>0');
else if k(i)==0.0001                 %若K=0.0001時，系統不穩定，判定系統不穩定
      disp('系統在K>0取值范圍不穩定');
else if k(i)==1.0001                 %若K=0.0001時，系統不穩定，在[0.001，1.001]
   k=0.0001:0.001:k(i);nn=length(k);            %區間搜索使系統臨界穩定的K值
   for i=1:nn
      p=chdihshu(Gsys,Hsys,k(i));  
      ii=find(real(p)>=0);n=length(ii);
      if(n>0),disp('系統穩定臨界值');k(i);break
   end;end
else                                      %若K>0.001時，系統才不穩定，在[K-1,K+1]
   k=k(i)-1:0.001:k(i)+1;nn=length(k);          %區間搜索使系統臨界穩定的K值 
   for i=1:nn
      p=chdihshu(Gsys,Hsys,k(i));    
      ii=find(real(p)>=0);n=length(ii);
      if(n>0),disp('系統穩定臨界值');k(i);break
   end,end
end;end;end
K=k(i);                                              %返回系統穩定臨界值K