帮我写一个受迫振动的matlab代码

以下是一个用 Matlab 编写的简单的受迫振动模拟程序：

%% 受迫振动模拟程序

% 物理参数
m = 1;     % 质量
k = 1;     % 弹性系数
w = 2*pi;  % 这里设定驱动频率为 1 Hz

% 初始状态
x0 = 0;   % 初位置
v0 = 0;   % 初速度

% 模拟参数设定
dt = 0.01;               % 时间步长
T = 20;                  % 总模拟时间
t = 0:dt:T;              % 时间序列
N = length(t);           % 时间点数
x = zeros(1,N);          % 位置序列
v = zeros(1,N);          % 速度序列
x(1) = x0; v(1) = v0;    % 初始化位置和速度

% 数值计算
for i = 2:N
    f = cos(w*t(i));  % 驱动力
    x(i) = x(i-1) + v(i-1)*dt;
    v(i) = v(i-1) - (k/m)*x(i-1)*dt + (1/m)*f*dt;
end

% 结果绘图
figure
plot(t,x)
title('Matlab 受迫振动模拟')
xlabel('Time [s]')
ylabel('y')

该程序利用欧拉法数值求解微分方程，模拟简谐运动的振幅随时间的变化。运行程序后，可得到如下图所示的振动曲线：