你这个并不是非线性系统

非常抱歉，我误解了你的问题。以下是一段用Newmark-beta积分法求三自由度非线性弹簧振子系统在随机激励下的响应的MATLAB代码：

function [t,y] = NonlinearSpringOscillator(beta,gamma,m,k,dt,tmax,F)
% beta: Newmark-beta方法的参数
% gamma: Newmark-beta方法的参数
% m: 质量矩阵
% k: 刚度矩阵
% dt: 时间步长
% tmax: 最大时间
% F: 随机激励力

n = tmax/dt; % 时间步数
t = linspace(0,tmax,n+1)'; % 时间向量
y = zeros(n+1,6); % 位移和速度矩阵
y(1,:) = [0;0;0;0;0;0]; % 初始条件

for i = 1:n
    % 计算加速度向量
    a = (F(i,:)'-k*y(i,1:3)'.^3)/m;
    
    % 计算位移和速度向量
    if i == 1
        y(i+1,1:3) = y(i,1:3) + y(i,4:6)*dt + 0.5*a*dt^2; % 第一步采用欧拉法
        y(i+1,4:6) = y(i,4:6) + a*dt;
    else
        y(i+1,1:3) = (m+beta*dt^2*k)\(F(i,:)'+k*(y(i,1:3)'+dt*y(i,4:6)'+0.5*dt^2*((1-2*beta)*y(i,1:3)'+2*beta*y(i-1,1:3)'.^3))); % 计算位移
        y(i+1,4:6) = gamma*dt*(y(i+1,1:3)-y(i,1:3))/dt + (1-gamma)*y(i,4:6) + dt*(1-gamma)*a'; % 计算速度
    end
end

end

这里假设非线性弹簧振子系统的刚度矩阵是按照 $k_{ij} = k_0 (x_{ij})^2$ 的形式计算的，其中 $x_{ij}$ 表示第 $i$ 个质点和第 $j$ 个质点之间的位移变量，$k_0$ 是弹簧的初始刚度系数。在计算加速度向量时，我们采用了 $k_{ij} = k_0 (x_{ij})^3$ 的形式，使得系统成为非线性的。