Vectorizing solve() for solving dynamical system for different parameter values

I will try to quickly explain the problem in brief and then pose the problem.

% Creating symbolic variables
% x0 y0 omega0 E0 are physical quantities
% Pm0 is a parameter generally in the range [0 1]
syms x0 y0 omega0 E0 Pm0
assume([x0 y0 omega0 E0], 'real');
assume(Pm0 ~= 0);
% x0 and y0 are related by
assumeAlso(x0*x0 + y0*y0 == 1);

% Equations for the dynamical system, (equating to 0 to find equilibrium
% points)
eqn1 = -y0*omega0 + x0*(1-x0*x0-y0*y0) == 0;
eqn2 = x0*omega0 + y0*(1-x0*x0-y0*y0) == 0;
eqn3 = -3.33*E0*y0 - 0.66*omega0 + 3.33 * Pm0 == 0;
eqn4 = 0.5*x0 - E0 + 0.5 == 0;
 
% solving the equations
eqn_final = [eqn1, eqn2, eqn3, eqn4];
variable_final = [x0, y0, omega0, E0];

a = solve(eqn_final, variable_final, "ReturnConditions", true);

% PmSol is a vector that contains different values for the parameter Pm0
PmSol = linspace(0.2, .5, 100);
% For each value of Pm0 from PmSol we need to solve the equations and find
% equilibrium pts.

% Then we need store the solutions for x0 to E0 for each Pm value.
% create vector of Fixed Point locations for each Pm val in PmSol
Fixed_points_locations = ones(2,length(PmSol)*4)*0;
% For each Pm value we will get either 2 solutions or no solutions
% if we get 2 solutions then we will store a vector like:
% [x0_1, y0_1 omega0_1 E0_1; x0_1, y0_1 omega0_1 E0_1] in above matrix. If
% no solutions then we will store zeros(2 by 4) in above matrix.

The PROBLEM:

% THE PROBLEM:
% THE AIM: the goal is to make the steps from here on out as fast as
% possible. This will become a part of a bigger problem needs to be
% lightning fast!


a

a = struct with fields:
            x0: [1×1 sym]
            y0: [1×1 sym]
        omega0: [1×1 sym]
            E0: [1×1 sym]
    parameters: [1×1 sym]
    conditions: [1×1 sym]

a.parameters

ans = 

x

a.conditions

ans = 

Pm_val = PmSol;
% substitute all Pm values at once in a.conditions to get equation for
% parameter
subs_new = subs(a.conditions, Pm0, Pm_val);  
subs_new(1,1)

ans = 

% Now all we need to do is solve each equation in subs_new (each of which
% looks like the above) get the parameter values, combine that with Pm_val
% and substitue in the expressions for x0 to E0 to get equilibrium points
% for all values of Pm

% THE ISSUE: HOW TO VECTORIZE THIS? (the following for loop)

for j = 1:length(PmSol)
    % solve the each equation in subs_new for the parameter
    parameter_val = solve(subs_new(j));
    % make a vector of same dimensions as parameter vals with value
    % Pm_val(j)
    Pm_new = parameter_val * 0 + Pm_val(j);
    
    % substitute parameter_val and Pm_new to get Equilibrium points
    T1 = {[parameter_val], [Pm_new]};
    T2 = {a.parameters, Pm0};
    xFP = subs(a.x0, T2, T1);
    yFP = subs(a.y0, T2, T1);
    omegaFP = subs(a.omega0, T2, T1);
    EFP = subs(a.E0, T2, T1);
    
    [rows, cols] = size([xFP, omegaFP, EFP]);
    % if no solutions add a zero matrix for the equilibrium point locations
    if rows == 0
        Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = ones(2, 3)*0;
        continue;
    end
    Fixed_points_locations(1:2,3*(j-1)+1:3*(j-1)+3) = [xFP, omegaFP, EFP];

end