Amplitude vs frequency curve

Hey guys ! I would like your help to plot this frequency response function.

Just so you guys understand, I'm trying to plot the frequency response curve through ode45. I have a one degree of freedom system with its constants and I'm driving this system with a sinusoidal force that will do a frequency sweep.

I am sending the analytic result, as for my code attempt. The curve using the fft appears to be coherent, but the y axis has a very different unit.

function simulation
clear all
close all
clc

m = 0.086; % kg
k = 166.3629; % N/m
c = 0.1664 ; % N.s/m
F =  0.6385;% N

% Time
Fs=400;
tspan =0:1/Fs:250-1/Fs;
f0 = 4; % Hz
f1 = 9; % Hz

% analytical solution
w = f0:0.2:f1; %Hz 
d1 = (k - m*(w*2*pi).^2).^2 + (c*w*2*pi).^2;
X = F./sqrt(d1);

% IC

x0 = 0; v0 = 0;
IC2 = [x0;v0];

% numerical integration
[time2,state_values2] = ode45(@h,tspan,IC2);
x = state_values2(:,1);

figure(1)
plot(time2,x),xlabel('time(s)'),ylabel('displacement(m)')

figure(2)
n=ceil(log2(length(x))); 
fx=fft(x,2^n);
fx=2*fx/length(x); % This operation is Adjusting the Magnitudes
f=(Fs/2^n)*(0:2^(n-1)-1);
plot(f,abs(fx(1:2^(n-1))),w,X,'-.k'), xlabel(' Frequency (Hz)'), ylabel(' X (m)');
l1 = ' using fft';
l2 = ' analytical solution';
legend(l1,l2);
xlim([f0 f1])
end


function sdot1 = h(t,x)

m = 0.086; % kg
k = 166.3629; % N/m
c = 0.1664 ; % N.s/m

f0 = 4; % Hz
f1 = 9; % Hz
F =  0.6385;% N
a = (f1 - f0)/250;
sdot1 = [x(2);
    (F.*sin(2*pi*(a*t/2 + f0)*t) - c.*x(2) - k*x(1))/m];

end

[Scott MacKenzie]

Is there a question herer?

You might consider providing more details, such as examples of the current output and code that can be executed to generate such output.

[José Anchieta de Jesus Filho]

First, thank you for your interest in helping me. Well, what my code should do is find the displacement amplitude as a function of each frequency, however, my fft seems to have the wrong units and I can't see what I'm doing wrong.

As you can see in the screenshot, when running the code, the fft is not in the same unit on the y-axis as the analytic solution. I would like some help to get this code to run correctly.

[Scott MacKenzie]

Is the y-axis magnitude really that important? Both solutions identify the signal at 7 Hz. Isn't that the key result? Sorry, if I'm completely off base here; this view is just based on my limited experience with spectrum analyses. I also notice that you resized the magnitudes after applying the fft. If you do so before, the result is closer to what you are looking for. Good luck.

x = rescale(x,-1, 1);
fx=fft(x,2^n);

[José Anchieta de Jesus Filho]

Thanks again for your help, but the y-axis is pretty important. I'm actually trying to plot the amplitude as a function of the excitation frequency. The analytical solution represents the amplitude as a function of the excitation frequency that this system will have, in steady state.