Running a function in “for” cycle, containing vpasolve()

I have some function in my script that solves a system of nonlinear algebraic equations

function Lam = lambdas1(x1, x2, x3, x4, w1, w2, g, tau, eta, zeta, sigma, eps, kap, gam_1, gam_2)

lam = 0.5 * atanh(-2*g / (w1 + w2 - 2*zeta - 2*eta));
mu = cosh(lam); nu = sinh(lam);
Om = -2*zeta - 2*eta*(nu/mu)^2 + 2*g*(nu/mu) + w2 * (nu/mu)^2  + w1;

A1 = mu^2*w1 + nu^2*w2; A2 = nu^2*w1 + mu^2*w2;
B = mu*nu*(w1 + w2); C = g*(mu^2 + nu^2) - 2*mu*nu*(zeta + eta);
F1 = 2*(zeta*mu^2 + eta*nu^2 - g*mu*nu);
F2 = 2*(zeta*nu^2 + eta*mu^2 - g*mu*nu);

Lam = vpasolve(-kap*(x1)^3 - kap*x2*(x1)^2 + 2*mu*tau*Om*(x1)^2 - 0.5*gam_1*x1 + ...
    (A1-F1)*x2 + (B-C)*x4 + 0.5*mu*tau*Om == 0, ...
    -(kap*(x1)^2*x2 + kap*(x2)^3 + 8*Om*tau^2*(x1)^3 + 4*mu*Om*tau*x1*x2 + ...
    (A1+F1+6*tau^2*Om)*x1 + 0.5*gam_1*x2 - (B+C)*x3 - nu*sigma + eps*mu) == 0, ...
    (B+C)*x2 - 0.5*gam_2*x3 + (A2+F2)*x4 == 0, ...
    -((C-B)*x1 + (A2-F2)*x3 + 0.5*gam_2*x4 + mu*sigma - eps*nu) == 0, ...
    [x1, x2, x3, x4]);

end

Where are real symbolic numbers which I define in the script (outside the sunction). I need to run this function in the cycle like this one (I have a vector ov values of zeta and I wish to solve a system for each of them)

syms x1 x2 x3 x4 real  

for k = 1 : some_number   

lam = 0.5 * atanh(-2*g / (w1 + w2 - 2*zeta(k) - 2*eta));
mu = cosh(lam); 
nu = sinh(lam);
Om = -2*zeta(k) - 2*eta*(nu/mu)^2 + 2*g*(nu/mu) + w2 * (nu/mu)^2  + w1;

A1 = mu^2*w1 + nu^2*w2; A2 = nu^2*w1 + mu^2*w2;
B = mu*nu*(w1 + w2); C = g*(mu^2 + nu^2) - 2*mu*nu*(zeta(k) + eta);
F1 = 2*(zeta(k)*mu^2 + eta*nu^2 - g*mu*nu);
F2 = 2*(zeta(k)*nu^2 + eta*mu^2 - g*mu*nu);

My_roots = vpasolve(-kap*(x1)^3 - kap*x2*(x1)^2 + 2*mu*tau*Om*(x1)^2 - 0.5*gam_1*x1 + ...
    (A1-F1)*x2 + (B-C)*x4 + 0.5*mu*tau*Om == 0, ...
    -(kap*(x1)^2*x2 + kap*(x2)^3 + 8*Om*tau^2*(x1)^3 + 4*mu*Om*tau*x1*x2 + ...
    (A1+F1+6*tau^2*Om)*x1 + 0.5*gam_1*x2 - (B+C)*x3 - nu*sigma + eps*mu) == 0, ...
    (B+C)*x2 - 0.5*gam_2*x3 + (A2+F2)*x4 == 0, ...
    -((C-B)*x1 + (A2-F2)*x3 + 0.5*gam_2*x4 + mu*sigma - eps*nu) == 0, ...
    [x1, x2, x3, x4]);

Solution_I_need(k) = My_roots.x1(1) + 1i*My_roots.x2(1);
   
end

But even for two iterations, it takes an extra-long time. Could you tell me please, how to deal with it properly?

[KSSV]

Why you want to solve for every iteration of loop? You can solve it for once, you end up with an expression; in this expression you subtstitute your values.

[Bogdan MP]

I'm solving numerically, this function would give me a set of roots for given parameters.

[Walter Roberson]

I think we will need some values to test with.

[Bogdan MP]

For example,

w1 = 2; 
w2 = 4;
tau = 1;

zeta = 0.1;
eta = 0.1;

sigma = 0;
eps = 0;

g = 0.3;

kap = 0.1;
gam_1 = 0.1; 
gam_2 = 0.1;

It runs easily and fast out of the cusle, but something ruins when I put it inside

[Walter Roberson]

I would suggest to you that it is not the fact that it is in a loop, but rather that your loop takes you to other zeta values that are more difficult to solve.

[Bogdan MP]

I've already thought so. But you may simply consider, let's say,

zeta(1) = 0.1;
zeta(2) = 0.2;

and calculate it out of a loop and inside. The first way is fast, the second is endless.

[Bogdan MP]

Yes, it looks like it's a problem with the system itself.

[Walter Roberson]

I am not clear at the moment why there is a problem. If you try 0.4 then it is very slow immediately. But when I use symbolic calculations, the forms I get with general zeta values can be solved fairly quickly numerically for any given zeta value, as it turns out to be just a degree 9 polynomial.

Speaking of degree 9 polynomials: it might be better if you supplied initial values for the vpasolve()

[Bogdan MP]

Did you simply run putting zeta as a symbolic variable? Because when I was trying to obtain the general result this way it was too long and I aborted running.