How do I plot a graph ode45

Question

0 votes

%%
clc
clear all
A1 = 0.015; %Excitation Amplitude
A2 = 0.035;
A3 = 0.06;
l_r = 0.617; %Wave length of the road
v = 0.04107647433; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
k_l = 26400; %Linear stiffness
m = 483; %Mass
d = -0.473; %Stretching condition
l = 0.946; %Length of the spring
k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
fig = figure();
ax = axes();
hold(ax); 
% view([-53 33]);
grid on
A_array = linspace(0,5,100); %Excitation Amplitude
T = 150;
x0 = [0,0];
xval = zeros([],1) ; 
yval = zeros([],1) ; 
for i=1:numel(A_array)
    A = A_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,[100,T],x0);
    A = A1*sin(Om*t)+A2*sin(2*Om*t)+A3*sin(3*Om*t);
    Response_amp = (max(x(:,1)) - min(x(:,1)))/2;
%    xval(i) = Om/(2*pi) ; 
    xval(i) = A; 
    yval(i) = Response_amp ; 
    
    % plot(Om, Response_amp, '.');
end
plot(xval,yval) ; 
xlabel('A (m)')
ylabel('Response Amplitude (m)')
ax = gca;
ax.XAxis.FontSize =16;
ax.YAxis.FontSize =16;