How to remove the time dependence of an equation?

I am usign the below code for plotting the final equation, eq_14 using the fimplicit function but the final equation has the form "eq_14(t)" and it has variables Omega and A as well. This is causing problem and hence I want the change "eq_14(t)=>eq_14". Please help me with this..

syms y(t) z0 xi delta alpha Omega ohm phi A theta p1 p2 p3 p4 p5 p6 p7 p8 

%%%%%%%%%% 5th order Equation %%%%%%%%%%%%%%
eq_1 = diff(y,t,t) + 2*xi*diff(y,t) + ...
    (p1*y^5 + p2*y^4 + p3*y^3 + p4*y^2 + p5*y^1 + p6*y^0 + p7) - z0

eq_1(t) = 

%%%%%%%%%% 5th order Equation %%%%%%%%%%%%%%
% eq_1 = diff(y,t,t) + 2*xi*diff(y,t) + ...
%     (p1*y^7 + p2*y^6 + p3*y^5 + p4*y^4 + p5*y^3 + p6*y^2 + p7*y + p8) - z0


%% === Substitution of assumed harmonic motion ===
xa = A*cos(Omega*t + phi)        % assumed displacement

xa = 

eq_2 = subs(eq_1, y, xa)        % substitute y = A*cos(Omega*t + phi)

eq_2(t) = 

eq_3 = subs(eq_2, (Omega*t + phi), theta)  % replace (Omega*t + phi) -> θ

eq_3(t) = 

%% === Expand and simplify ===
eq_4 = expand(eq_3)

eq_4(t) = 

eq_5 = combine(eq_4, 'sincos')   % neatly group sin/cos terms

eq_5(t) = 

disp('Combined equation in sin(theta) and cos(theta):')

Combined equation in sin(theta) and cos(theta):

% pretty(eq_5)

%% === Extract only coefficients of sin(theta) and cos(theta) ===
% --- sin(theta) coefficient ---
eq_6_sin = simplify(subs(eq_5, sin(theta), 1) - subs(eq_5, sin(theta), 0))

eq_6_sin(t) = 

%%%%%% === Display results ===
% --- cos(theta) coefficient ---
eq_7_cos = simplify(subs(eq_5, cos(theta), 1) - subs(eq_5, cos(theta), 0))

eq_7_cos(t) = 

disp('Coefficient of sin(theta):')

Coefficient of sin(theta):

% Collect both for clarity
eq_8 = simplify(collect(eq_6_sin, [A Omega]))

eq_8(t) = 

disp('Coefficient of cos(theta):')

Coefficient of cos(theta):

eq_9 = simplify(collect(eq_7_cos, [A Omega]))

eq_9(t) = 

% pretty(eq_sin)
% pretty(eq_cos)
% further sqauring and adding than will get FRF Quadrating equation 

% squaring eq_8 and eq_9 and z_0 and adding 
% Note - Z0^2 add in eq_10 because in motion of equion right side having
% exciation z0 i.e...(% 'Equation of motion: x'' + 2*xi*x' + f_QZS = -z0*cos(Omega*t)
eq_10=eq_8^2+eq_9^2-z0^2

eq_10(t) = 

eq_11 = expand(eq_10)

eq_11(t) = 

eq_12=simplify(eq_11)

eq_12(t) = 

eq_13=collect(eq_12,Omega)

eq_13(t) = 

eq_13=collect(eq_12,Omega)

eq_13(t) = 

eq_13=expand(eq_13/A^2)

eq_13(t) = 

eq_13=collect(eq_13, Omega)

eq_13(t) = 

Vxi=0.1;
Vz0=0.075;
Vp1= 0.0232;
Vp3= 2.5523e-4;
Vp5= 2.1850e-7;

eq_14 = subs(eq_13,[p1,p3,p5,xi,z0],[Vp1,Vp3,Vp5,Vxi,Vz0])

eq_14(t) = 

figure; hold on;

fimplicit(eq_14==0,"r","LineWidth",2);hold on; %,[0 5 0 0.8],

Error using fimplicit>singleFimplicit (line 185)
Input must be a function or functions of a single variable.

Error in fimplicit>@(f)singleFimplicit(cax,f,limits,extraOpts,args) (line 154)
hObj = cellfun(@(f) singleFimplicit(cax,f,limits,extraOpts,args),fn,'UniformOutput',false);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fimplicit>vectorizeFimplicit (line 154)
hObj = cellfun(@(f) singleFimplicit(cax,f,limits,extraOpts,args),fn,'UniformOutput',false);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fimplicit (line 128)
hObj = vectorizeFimplicit(cax,fn,limits,extraOpts,args);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

xlabel('Frequency ratio');                            % enter x label
ylabel('Amplitude'); 


%%%%% although I am able to plot this equation which is same as above but
%%%%% without any time dependance 

    eq13=Omega^4+...
        Omega^2*(-5/4*p1*A^4-3/2*p3*A^2+4*xi^2-2*p5) ...
        +9/16*p3^2*A^4 ...
        +25/64*p1^2*A^8 ...
        +p5^2 ...
        -z0.^2/A^2 ...
        +3/2*p3*p5*A^2 ...
        +5/4*p1*p5*A^4 ...
        +15/16*p1*p3*A^6;


    var = subs(eq13,[p1,p3,p5,xi,z0],[Vp1,Vp3,Vp5,Vxi,Vz0]);
    var = subs(var,[A,Omega],[y,x]);
    
    figure; hold on;
    fimplicit(var==0,"r","LineWidth",2);hold on;
    xlabel('Frequency ratio');                  
    ylabel('Amplitude'); 

Please suggest me a way to findout the solution for this.

[Walter Roberson]

syms y(t) z0 xi delta alpha Omega ohm phi A theta p1 p2 p3 p4 p5 p6 p7 p8 

%%%%%%%%%% 5th order Equation %%%%%%%%%%%%%%
eq_1 = diff(y,t,t) + 2*xi*diff(y,t) + ...
    (p1*y^5 + p2*y^4 + p3*y^3 + p4*y^2 + p5*y^1 + p6*y^0 + p7) - z0

eq_1(t) = 

%%%%%%%%%% 5th order Equation %%%%%%%%%%%%%%
% eq_1 = diff(y,t,t) + 2*xi*diff(y,t) + ...
%     (p1*y^7 + p2*y^6 + p3*y^5 + p4*y^4 + p5*y^3 + p6*y^2 + p7*y + p8) - z0


%% === Substitution of assumed harmonic motion ===
xa = A*cos(Omega*t + phi)        % assumed displacement

xa = 

eq_2 = subs(eq_1, y, xa)        % substitute y = A*cos(Omega*t + phi)

eq_2(t) = 

eq_3 = subs(eq_2, (Omega*t + phi), theta)  % replace (Omega*t + phi) -> θ

eq_3(t) = 

%% === Expand and simplify ===
eq_4 = expand(eq_3)

eq_4(t) = 

eq_5 = combine(eq_4, 'sincos')   % neatly group sin/cos terms

eq_5(t) = 

disp('Combined equation in sin(theta) and cos(theta):')

Combined equation in sin(theta) and cos(theta):

% pretty(eq_5)

%% === Extract only coefficients of sin(theta) and cos(theta) ===
% --- sin(theta) coefficient ---
eq_6_sin = simplify(subs(eq_5, sin(theta), 1) - subs(eq_5, sin(theta), 0))

eq_6_sin(t) = 

%%%%%% === Display results ===
% --- cos(theta) coefficient ---
eq_7_cos = simplify(subs(eq_5, cos(theta), 1) - subs(eq_5, cos(theta), 0))

eq_7_cos(t) = 

disp('Coefficient of sin(theta):')

Coefficient of sin(theta):

% Collect both for clarity
eq_8 = simplify(collect(eq_6_sin, [A Omega]))

eq_8(t) = 

disp('Coefficient of cos(theta):')

Coefficient of cos(theta):

eq_9 = simplify(collect(eq_7_cos, [A Omega]))

eq_9(t) = 

% pretty(eq_sin)
% pretty(eq_cos)
% further sqauring and adding than will get FRF Quadrating equation 

% squaring eq_8 and eq_9 and z_0 and adding 
% Note - Z0^2 add in eq_10 because in motion of equion right side having
% exciation z0 i.e...(% 'Equation of motion: x'' + 2*xi*x' + f_QZS = -z0*cos(Omega*t)
eq_10=eq_8^2+eq_9^2-z0^2

eq_10(t) = 

eq_11 = expand(eq_10)

eq_11(t) = 

eq_12=simplify(eq_11)

eq_12(t) = 

eq_13=collect(eq_12,Omega)

eq_13(t) = 

eq_13=collect(eq_12,Omega)

eq_13(t) = 

eq_13=expand(eq_13/A^2)

eq_13(t) = 

eq_13=collect(eq_13, Omega)

eq_13(t) = 

Vxi=0.1;
Vz0=0.075;
Vp1= 0.0232;
Vp3= 2.5523e-4;
Vp5= 2.1850e-7;

eq_14 = subs(eq_13,[p1,p3,p5,xi,z0],[Vp1,Vp3,Vp5,Vxi,Vz0])

eq_14(t) = 

figure; hold on;
fimplicit(eq_14(t)==0,"r","LineWidth",2);hold on; %,[0 5 0 0.8],
xlabel('Frequency ratio');                            % enter x label
ylabel('Amplitude'); 

%%%%% although I am able to plot this equation which is same as above but
%%%%% without any time dependance 

    eq13=Omega^4+...
        Omega^2*(-5/4*p1*A^4-3/2*p3*A^2+4*xi^2-2*p5) ...
        +9/16*p3^2*A^4 ...
        +25/64*p1^2*A^8 ...
        +p5^2 ...
        -z0.^2/A^2 ...
        +3/2*p3*p5*A^2 ...
        +5/4*p1*p5*A^4 ...
        +15/16*p1*p3*A^6;


    var = subs(eq13,[p1,p3,p5,xi,z0],[Vp1,Vp3,Vp5,Vxi,Vz0]);
    
    %at this point in your code, x is undefined, so guess about the
    %definition
    syms x(t)
    var = subs(var,[A,Omega],[y,x])

var = 

    
    figure; hold on;
    fimplicit(var==0,"r","LineWidth",2);hold on;
    xlabel('Frequency ratio');                  
    ylabel('Amplitude'); 

●

Warning: Error in state of SceneNode.
Unable to convert symbolic expression to double array because it contains symbolic function that does not evaluate to number. Input expression must evaluate to number.

fimplicit fails because var contains x(t) and y(t), both of which are undefined functions.

Answer 1

[Torsten]

Use

fimplicit(eq_14(t)==0,"r","LineWidth",2)

instead of

fimplicit(eq_14==0,"r","LineWidth",2)

Why do you use a time-dependent y(t) at all ? You never explicitly solve the differential equation in your code.

[Ron]

Firstly I want to thank you for this, I am doing so because I want the equation in the form of omega. If you have any other way of doing so then please let me know.

[Torsten]

solve(eq_13(t)==0,Omega)

Or substitute Omega^2 = new_variable.

This gives a quadratic equation in "new_variable", but should lead to the same result.