利用龙格-库塔法求解非线性振动方程的响应。给定非线性振动方程为mx 的二次导+ cx的一次导 + kx + k ₃x ³= f cos ( ωt )式中，m = 1；c = 0.04；k = 1；k₃ = 0.2；f = 3；ω= 0.9。求此MATLAB程序代码

龙格-库塔法（Runge-Kutta method），是一种数值积分方法，常用于解决微分方程组。对于非线性振动方程，我们可以将其转化为一阶常微分方程系统，然后应用四阶龙格-库塔（RK4）算法。以下是一个基本的MATLAB代码示例，它假设`x`和`dxdt`分别为状态变量和其导数：

function dxdt = nonLinearVibration(x, t)
% 非线性振动方程的定义
m = 1;
c = 0.04;
k = 1;
k3 = 0.2;
f = 3;
omega = 0.9;

% 计算右侧函数
rhs = [x(2); -m*x(1) - c*x(2) - k*x(1) - k3*x(1)^3 + f*cos(omega*t)];

% 返回导数向量
dxdt = rhs;
end

% 设置初始条件和时间范围
x0 = [0; 0]; % 初始位置和速度
tspan = [0 10]; % 时间范围从0到10秒

% 使用四阶龙格-库塔法求解
[t, x] = ode45(@nonLinearVibration, tspan, x0);

% 结果存储在矩阵x中，每一行对应一个时间点的解
% 可以通过 plot(t, x(:,1), 'b') 和 plot(t, x(:,2), 'r') 来分别绘制位移和速度随时间的变化

运行这段代码后，你需要查看`x`矩阵来获取各个时刻的状态。注意，实际使用时可能需要调整时间步长和边界条件，根据精度需求优化算法。