clc,clear all
A1 = 0.015;
A2 = 0.035;
A3 = 0.06;
l_r = 0.617;
v = 10;
P = l_r/v;
Om = 1/P*2*pi;
k_l = 26400;
m = 483;
d = -0.5;
l = 0.5;
k_s = -(k_l*(l-d))/(4*d);
f_n = sqrt(k_l/m)/(2*pi);
fig = figure();
ax = axes();
hold(ax);
view([45 85]);
grid on
l_array = 0.1:0.1:1;
T = 80;
x0 = [0,0];
for i=1:numel(l_array)
l = l_array(i);
f = @(t,x) [ x(2); ...
-(2*k_s*(x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))* ...
(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2) - l))/ ...
(m*(sqrt((l-d)^2 + (x(1)-(A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t)))^2))) ];
[t, x] = ode45(f,[0,T],x0);
plot3(t, l*ones(size(t)), x(:,1), 'LineWidth', 1);
end
xlabel('Time (s)', 'Fontsize', 14)
ylabel('Length (m)', 'Fontsize', 14)
zlabel('Response Amplitude (m)', 'Fontsize', 14)
ax = gca;
ax.XAxis.FontSize =14;
ax.YAxis.FontSize =14;
ax.ZAxis.FontSize =14;