Hi @Jerry Yeung
Basically, there is nothing wrong with your original code. The reason that the spring response is not affected by the Driving Force is because the magnitude of the Driving Force (5 N) to too small (almost negligible) relative to the Spring's initial Restoring Force 
, when 
is a very large value and
.
, when is a very large value and.So, either you select a smaller initial value for 
where 5 N is relatively significant, or you design a higher magnitude for the driving force. Two cases are included in the code. The Driving Force in Case #2 is state-dependent so that the initial values can be arbitrarily selected.
where 5 N is relatively significant, or you design a higher magnitude for the driving force. Two cases are included in the code. The Driving Force in Case #2 is state-dependent so that the initial values can be arbitrarily selected.On practical implementation, you need to check if the initial position is realistically selected at 1 meter. If it is true, then you also need to check if you can find a machine that delivers that amount of large driving force.
is realistically selected at 1 meter. If it is true, then you also need to check if you can find a machine that delivers that amount of large driving force.tend = 3e-13; % end time
rtol = 3e-14; % relative tolerance
atol = 1e-28; % absolute tolerance
options = odeset('RelTol', rtol, 'AbsTol', atol);
% p0 = [1.6e-30 0]; % initial position and velocity CASE #1
p0 = [1 0]; % initial position and velocity CASE #2
[t, p] = ode45(@(t, y) SpringFunction(t, y), [0 tend], p0, options);
plot(t, p(:, 1)), % plot the first column (x) vs. time
grid on, xlabel('t')
function dydt = SpringFunction(t, y)
c0 = 299792458; % speed of light
lase_lambda = 1064e-9; % wavelength
lase_freq = c0/lase_lambda;
w0 = 2*pi*lase_freq; % Oscillation frequency
tau2 = 500e-6;
c = 1/tau2*1e10; % Damping coefficient
% f = 5*cos(t*w0); % Driving force CASE #1
wc = 6e13; % design parameter for Case #2
f = (c - 2*wc)*y(2) + (w0^2 - wc^2)*y(1); % CASE #2
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = f - c*y(2) - (w0^2)*y(1);
end
