在MATLAB中写出图中的二阶微分方程，通过四阶龙格库塔法求解x(t),y(t)的离散解，离散解结果无误，但是根据离散解通过写的傅里叶变换程序求功率谱如左图，始终无法画出右图中的结果。而且画出来的结果一直受计算时长的影响较大，但参考文献中没有写出。想请问应该写怎么样的程序才能得出右图的功率谱。
clear;
clc
fwind=0.01144;
R1=3.17;
R3=-1.61;
R5=-2.76;
g=9.81;
I=5.37*10^10;
delta=1.47*10^8;
r1=delta*g*R1/I;
r3=delta*g*R3/I;
r5=delta*g*R5/I;
%options=odeset('maxstep',0.01);%%时间步长调整
%[t,y]=ode45(@vdp1,[-30 30],[0.875;0],options);%%四阶龙哥库塔法求解微分方程组，y(:,1)即x的解，y(:,2)即y的解
odefun2 = @(t,y)[y(2);-(r1*y(1)+r3*y(1)^3+r5*y(1)^5)+fwind*(1-y(1)^2/2+y(1)^4/24)]; % 化二阶常微分方程为二阶常微分方程组
ts = -30; %计算开始时间
te = 30; %计算结束时间
h = 0.001; %时间步长
t = ts:h:te; %计算数据点
y = zeros(length(t),2);
y(1,:)=[0.875,0]; %参数初值
for n = 1:length(t)-1
k1 = odefun2(t(n),y(n,:)');
k2 = odefun2(t(n)+h/2,y(n,:)'+h/2*k1); %这里要特别注意同型向量相加
k3 = odefun2(t(n)+h/2,y(n,:)'+h/2*k2);
k4 = odefun2(t(n)+h,y(n,:)'+h*k3);
y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';%四级四阶Runge-Kutta公式
end
%%傅里叶变换
N=length(y(:,2));
y_OMEGA=fft(y(:,2));
f=(0:N-1)/(N*h);%%采样频率
omega_0=(delta*g*R1/I)^0.5;
omega=f*2*pi;
OMEGA=omega/omega_0; %求出相应横坐标Ω
S_y=(abs(y_OMEGA)).^2./(N); %Sy(Ω)
plot(OMEGA,S_y);
axis([0 1.5, -inf,inf]);
