a=0.4;
p=0.2;
b=1.29;
w0=0.8;
T=2*pi/w0;
[t,y]=ode45('liangfun',[0:T/300:400*T],[0.1,0],[],a,b,p,w0);
%畫功率譜圖
%****************************************************
Y=fft(y(:,1));
Y(1)=[];
n=length(Y);
power=abs(Y(1:n/2)).^2/n^2;
freq=300*(1:n/2)./n;
power=10*log(power);
figure(1)
plot(freq,power,'k')
%axis([0,2,0,0.35])
xlabel('頻率/Hz');
ylabel('功率(10*log(power))/ω');
%*****************************************************
%畫poincare映射圖
%*****************************************************
figure(2)
plot(y(12000:300:end,1),y(12000:300:end,2),'k.');
xlabel('δ/rad'),ylabel('ω/(rad/s)');