Frequency vs Amplitude using ODE45 of system of equations?

I have a system of equations and I want to plot the frequency vs amplitude curves for the system.

this is my function

function xdot= multifunction(t,x)
 
global    alfa Meo v  eta zeta P vt psi
 
xdot = zeros(8,1);
 
 
xdot(1)=x(2); 
 
xdot(2)= - ((1/2)*pi^4 *(zeta+1)*x(1)) - Meo *x(2) - (1/8)* eta * (-x(6)*x(1)^2 + 4*x(6)*x(1)*x(5) + x(6)*x(3)*(x(3)+2*x(7)) ...
+ 2*x(8)*(-x(1)*x(3) + x(3)*x(5) + 2*x(1)*x(7)) + 2*x(4)*(x(1)*(2*x(3)-x(7))+ x(5)*(x(3)+x(7))) + x(2)*(3*x(1)^2 - 2*x(1)*x(5) ...
+ 2*(x(3)^2 + x(5)^2 - x(3)*x(7) + x(7)^2))) ...
- (1/8)*pi^4*x(1)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) ...
- (1/2)*alfa*v*pi*(-2*x(2)*x(3) + 2*x(6)*x(3) + x(4)*(x(1) - 3*x(5)) + 3*x(8)*x(5) - 4*x(6)*x(7)) + ((1/2)*pi^2 * x(1)*(alfa^2 * v^2 - P)) ...
- (psi*vt)*((4/15)*(5*x(3) + 2*x(7)));
 
 
 
xdot(3)=x(4); 
 
xdot(4)= - ((8)*pi^4 *(zeta+1)*x(3)) - Meo *x(4) - (1/8)* eta * (x(8)*(-x(1)^2 + 2*x(1)*x(5) + x(5)^2 + 4*x(3)*x(7)) ...
+2*(x(2)*x(1)*(2*x(3) - x(7)) + x(4)*x(1)*(x(3) + x(7)) + x(2)*x(5)*(x(3) + x(7))+ x(6)*x(5)*(2*x(3) + x(7))) ...
+ x(4)*(2*x(1)^2 + 2*x(1)*x(5) + 3*x(3)^2 + 2*x(5)^2 +  2*x(7)^2)) - (1/2)*pi^4*x(3)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) ...
- (1/2)*alfa*v*pi*( x(7)*x(1) + 2*x(8)*x(3) + x(2)*(x(1) - 3*x(5)) - 4*x(4)*x(7)) + (2*pi^2 * x(3)*(alfa^2 * v^2 - P)) ...
- (psi*vt)*((4/15)*(5*x(1) - 9*x(5)));
 
 
 
xdot(5)=x(6); 
 
xdot(6)= - ((81/2)*pi^4 *(zeta+1)*x(5)) - Meo *x(6) - (1/8)* eta * (2*x(6)*x(1)^2 + 3*x(6)*x(5)^2 + 2*x(8)*x(3)*(x(1)+x(5)) ...
+ 2*x(4)*x(3)*(x(1)+2*x(5)) + 4*x(8)*x(5)*x(7) +2*x(4)*x(7)*(x(1)+x(5)) + 2*x(6)*(x(3)^2 + x(3)*x(7)+x(7)^2) + x(2)*(-x(1)^2 + 4*x(1)*x(5) + x(3)*(x(3)+2*x(7)))) ...
- (9/8)*pi^4*x(5)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) - (1/2)*alfa*v*pi*( x(4)*x(1) + x(8)*x(1) + 2*x(2)*(x(3) - 2*x(7)));
+ ((9/2)*pi^2 * x(5)*(alfa^2 * v^2 - P)) - (psi*vt)*((12/35)*(7*x(3) - 10*x(7)));
 
 
 
xdot(7)=x(8); 
 
xdot(8)= - ((128)*pi^4 *(zeta+1)*x(7)) - Meo *x(8) - (1/8)* eta * ((2*(x(2)*x(3))*(-x(1) + x(5)) + x(6)*x(3)*(x(1)+x(5)) + x(8)*(x(1)^2 + x(3)^2 +x(5)^2)) ...
+ 4*x(7)*(x(2)*x(1) + x(6)*x(5)) + 3* x(8)*x(7)^2 + x(4)*(-x(1)^2 + 2*x(1)*x(5) + x(5)^2 + 4*x(3)*x(7))) ...
- 2*pi^4*x(7)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) - (1/2)*alfa*v*pi*(x(6)*x(1) + 2*x(4)*x(3) + 3*x(2)*x(5));
+ (8*pi^2 * x(7)*(alfa^2 * v^2 - P)) - (psi*vt)*((8/105)*(7*x(1) + 45*x(5)));
 
end 

all material properties and functions are specified and coded in another script. I only need the ODE45 code to execute and get the frequency vs amplitude curves

The script I wrote is to plot the Phase portrait and the wave form.

clc;clear all;close all;format long;
global alfa Meo v vt eta zeta P psi

E=1.1e12;
row=1.3e3;
G=280e9;
L=60e-9;
h=1e-9;
w= 10*h;
gap= 5*h;
I=(1/12)*w*(h^3);
Ar=w*L;
v= 2;
vt = 1;
Meo1=200e-8;
Po = 1e5;
Pa = 1e7;
lambda=64e-9;
lam = (lambda*Po)/(Pa*gap);
cpp = (Meo1*(L^2)*(w^3))/((gap^3)*sqrt(E*I*Ar*row));
Pr=100;
ssf = h;
eta1 = 2e7; 
r = sqrt(1/12)*h;
 
 
 
zeta = ((L^2)*G*ssf)/(E*(r^2));
alfa = r/L;
P = (Pr*L^2)/(E*(r^2));
eta = (eta1*(L^2)*sqrt(row/E)/(Ar*row));
Meo = (cpp*(L^2)*sqrt(row/E)/(Ar*row*r));
psi = L*sqrt((row*Ar)/(E*I));
 
%% ODE Solver
t0=0;
Tt=2*pi/81.8
Tf=500*Tt
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,x]=ode45(@multifunction,0:Tt/100:Tf,[.01 0 0 0 0 0 0 0],options);
n=length(t);
as=fix(0.98*n);
t=t(as:n);
x1=x(as:n,1);
x2=x(as:n,2);
 
%% Results
figure(1)
set(gcf,'Color',[1,1,1])
 
subplot(1,2,1);hold on; 
plot(t,x1,'k','linewidth',0.1)
grid on
title (' Wave Form','FontSize',12)
xlabel('Time','FontSize',12)
ylabel('$q1$', 'interpreter', 'latex','FontSize',20)

%%
subplot(1,2,2);hold on; 
plot(x1,x2,'k','linewidth',0.1)
grid on
title ('Phase Portrait','FontSize',12)
xlabel('$q1$', 'interpreter', 'latex','FontSize',20)
ylabel('$\dot q1$', 'interpreter', 'latex','FontSize',20)

Please help!