Faster processing time for code

I've got some code that operates how I want it to. It takes some input values and calculates the quartic roots of an equation containing those input values. The issue is that the code takes ~6hrs to fully run. I have shortened the range of values I'm testing, but it's still a hassle if I want to check results, seeing as the code takes so long to run. I was looking for any advice or alternative methods on running this code more quickly. Thanks for your assistance!

clear
clc

%%CONSTANTS%%
c = 2.998e8; %ms-1, speed of light
e = 1.602e-19; %C, electric charge
m_p = 1.673e-27; %kg, proton mass
m_e = 9.109e-31; %kg, electron mass
K = 1.381e-23; %JK-1, Boltzmann constant
eps_0 = 8.854e-12; %F/m, permittivity of free space


T_e = 1e5; %electron temperature in Kelvin
n_0 = 10e6; %density of 10 particles/cubic cm 

V_s = ((K*T_e)/m_p)^(1/2); %ion sound speed
v_b = 0.8 * V_s; %beam speed, m/s

w_pe = ((n_0 * e^2)/(eps_0 * m_e))^(1/2);
w_range = [1e-4:1e-4:1e-1];
w = w_range .* w_pe;

V_e = sqrt((K*T_e)/m_e);
lambda_D = w_pe/V_e;

klambda_D = [1e-4:1e-4:1e-1];

k = [klambda_D./lambda_D]

n_bn_0 = 1e-3;

%%FINISH CONSTANTS%%

A = (k.^2 .* V_s.^2)./(1 + klambda_D.^2);

a = 1;
b = -2 .* k .* v_b^2;
d = (k.^2 .* v_b.^2) - A - n_bn_0 .* A;
f = 2 .* A .* k .* v_b;
g = -A .* k.^2 .* v_b.^2;

N = numel(k);
R = nan(4,N);

for poly = 1:N % loop over indices
    poly = [a, b, d, f, g]
    R = roots(poly)
end

format long 

Roots = array2table(R)
R_real = real(R)
R_imag = imag(R)

[Dyuman Joshi]

You are just running the same loop over and over again.

Do you want to get roots corresponding to each value in k?

Answer 1

[VBBV]

did you try changing step size ?

clear
clc

tic
%%CONSTANTS%%
c = 2.998e8; %ms-1, speed of light
e = 1.602e-19; %C, electric charge
m_p = 1.673e-27; %kg, proton mass
m_e = 9.109e-31; %kg, electron mass
K = 1.381e-23; %JK-1, Boltzmann constant
eps_0 = 8.854e-12; %F/m, permittivity of free space


T_e = 1e5; %electron temperature in Kelvin
n_0 = 10e6; %density of 10 particles/cubic cm 

V_s = ((K*T_e)/m_p)^(1/2); %ion sound speed
v_b = 0.8 * V_s; %beam speed, m/s

w_pe = ((n_0 * e^2)/(eps_0 * m_e))^(1/2);
w_range = [1e-4:1e-4:1e-1];
w = w_range .* w_pe;

V_e = sqrt((K*T_e)/m_e);
lambda_D = w_pe/V_e;

% did you try increasing step size
klambda_D = [1e-4:1e-2:1e-1];

k = [klambda_D./lambda_D]

k = 1×10

0.0007 0.0697 0.1387 0.2078 0.2768 0.3458 0.4148 0.4839 0.5529 0.6219

n_bn_0 = 1e-3;

%%FINISH CONSTANTS%%

A = (k.^2 .* V_s.^2)./(1 + klambda_D.^2);

a = 1;
b = -2 .* k .* v_b^2;
d = (k.^2 .* v_b.^2) - A - n_bn_0 .* A;
f = 2 .* A .* k .* v_b;
g = -A .* k.^2 .* v_b.^2;

N = numel(k);
R = nan(4,N);

for poly = 1:N % loop over indices
    poly = [a, b, d, f, g];
    R = roots(poly);
end

toc

Elapsed time is 0.100909 seconds.

format long 

Roots = array2table(R);
R_real = real(R);
R_imag = imag(R);

[VBBV]

Note: You dont need a for loop in this case,

poly = [a, b, d, f, g];
R = roots(poly);

The above lines are enough. Also the fo loop seem to be used incorrectly