function [DX1,DX2,DP]=portraitzzy(x1,x2,h)
% PORTRAITZZY		采用"一步仿真"計算狀態變量斜率和狀態導數的二次方根
% x1,x2				分別給定"狀態平面"上的格點坐標
% h						給定積分計算采用的時間步長
% DX1,DX2				軌線斜率在狀態坐標軸上的投影長度
% DP					狀態導數向量的長度
opts=simset('solver','ode5','fixedstep',h);		%采用ode5定步長積分算法	<7>
n=length(x1);			
X1=zeros(n,n);X2=X1;									%預置空間
for ii=1:n;
   for jj=1:n;
      opts=simset(opts,'initialstate',[x1(ii),x2(jj)]);%設置狀態初值		<12>
      [t,x,y]=sim('exm09634_1',h,opts);	%步長為 h 的"一步仿真"		<13>
      dx1=x(2,1)-x1(ii);					%計算x1的變化率
      dx2=x(2,2)-x2(jj);					%計算x2的變化率
      L=sqrt(dx1^2+dx2^2);					%計算狀態軌線長度的變化率。
      Z(jj,ii)=L;
      if L>1.e-10								%若狀態軌線變化率大于"零"閾值
         DX1(jj,ii)=dx1/L;DX2(jj,ii)=dx2/L;	%計算各狀態變量的近似斜率	<19>
         										%注意下標次序。這是繪圖指令格式要求		
      end
   end
end
DP=Z/h;										%狀態導數向量的長度