Plotting within for loop

I've posted a couple of sections of code previously and gotten some really helpful advice! But I figure that if I put all my code in as I explain my issue this time maybe some intelligent soul can help me locate any issues I'll be having.

Essentially I have a cubic equation. I'm solving the cubic equation for given values of alpha_0, then separating each of the resulting roots into real and imaginary parts. The real and imaginary parts are then separately put into an equation and plotted against each other. This should result in one plot for each given value of alpha_0, with 6 lines for each of the real and imaginary root equations.

Currently whilst I can get the correct number of plots to appear, they are all identical and are not plotting the necessary values. I'd really appreciate it if anyone has any advice or corrections they could offer to correct my code? Thanks so much!

alpha_0 =  [-1,-0.5,-0.1,0,0.1,0.5,1];

loadconstants

k = w/V_s;

omega1 = (k .* V_s)./((1 + (klambda_D).^2).^(1/2))/w_pe;
omega2 = (k .* v_b)/w_pe;

N = numel(alpha_0);
R = nan(3,N);

for alpha = 1:N % loop over indices
    p = [1, 2*alpha_0(alpha), alpha_0(alpha)^2, -1/2];
    R(:,alpha) = roots(p);
end
Roots = array2table(R)
r = real(R)
i = imag(R)

rl = 1:numel(r);
im = 1:numel(i);

[numRows, numCols] = size(R)

hold on
grid on

lgd = legend('Location', 'northwest'); 
lgd.FontSize = 14;
lgd.Title.String = 'Disp relations';


for col = 1:numCols %for each column of R
    
    for r = r(:,1)
        delta_w_r = (r * (n_bn_0)^(1/3)) .* (k .* V_s)./((1 + (klambda_D).^2).^(1/2));
        omegar = ((k .* V_s)./((1 + (klambda_D).^2).^(1/2)) + delta_w_r)/w_pe;
        
        q = plot(omegar, klambda_D, '-');
        ylim([-0.1 0.1]);
        hold on
        grid on
        
        lgd = legend('Location', 'northwest'); 
        lgd.FontSize = 14;
        lgd.Title.String = 'Disp relations';
    end
    
    for i = i(:,1)
        delta_w_i = (i * (n_bn_0)^(1/3)) .* (k .* V_s)./((1 + (klambda_D).^2).^(1/2));
        omegai = ((k .* V_s)./((1 + (klambda_D).^2).^(1/2)) + delta_w_i)/w_pe;
        
        b = plot(omegai, klambda_D, '--');
        
    end
    h = figure;
    
    
    xlabel('k * \lambda_D')
    ylabel(' omega \omega')
    hold on
    
    
    
    saveas(h,sprintf('omega3%d.png',col));
    hold off
end    

Answer 1

[Voss]

for r = r(:,1)

This redefines the variable r. Before the loop executes, r is a 3-by-7 matrix. After the loop, r will be a 3-by-1 column vector and you won't be able to get any of the remaining columns that used to be in r, because r has been redefined to include only its first column.

Recall that a for loop iterates over the columns of what you give it, so if you say, "for index = valArray", then (from the documentation):

• valArray — Create a column vector, index, from subsequent columns of array valArray on each iteration. For example, on the first iteration, index = valArray(:,1). The loop executes a maximum of n times, where n is the number of columns of valArray, given by numel(valArray(1,:)). The input valArray can be of any MATLAB® data type, including a character vector, cell array, or struct.

Therefore, giving a for loop a column vector causes it to iterate only once, so the for loop is not even necessary. It's the same as saying r = r(:,1) and then executing the code inside the loop once. Looked at this way, of course you only get results corresponding to the first column of r.

The same problem applies to the other loop, which uses "for i = i(:,1)".

What you can do is get rid of those for statements and put the code inside those loops directly inside the outer loop over col, and that becomes the only loop you'll have. Something like this (I'm also renaming the variables r and i because for one thing it's not advisable to use i since it can be confused with the imaginary unit 1i):

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

[numRows, numCols] = size(R)

for col = 1:numCols %for each column of R
    
    % create a new figure
    h = figure;
    
    % calculate delta_w_r and omegar for this column of R_real, i.e., R_real(:,col)
    delta_w_r = (R_real(:,col) * (n_bn_0)^(1/3)) .* (k .* V_s)./((1 + (klambda_D).^2).^(1/2));
    omegar = ((k .* V_s)./((1 + (klambda_D).^2).^(1/2)) + delta_w_r)/w_pe;
    
    % plot
    q = plot(omegar, klambda_D, '-');
    ylim([-0.1 0.1]);
    hold on
    grid on
    
    % calculate delta_w_i and omegai for this column of R_imag, i.e., R_imag(:,col)
    delta_w_i = (R_imag(:,col) * (n_bn_0)^(1/3)) .* (k .* V_s)./((1 + (klambda_D).^2).^(1/2));
    omegai = ((k .* V_s)./((1 + (klambda_D).^2).^(1/2)) + delta_w_i)/w_pe;
    
    % plot
    b = plot(omegai, klambda_D, '--');
    
    % create a legend
    lgd = legend('Location', 'northwest'); 
    lgd.FontSize = 14;
    lgd.Title.String = 'Disp relations';
    
    % set axes labels
    xlabel('k * \lambda_D')
    ylabel(' omega \omega')
    
    % save the figure
    saveas(h,sprintf('omega3%d.png',col));
    
end

(Note that I also moved where you create the figures because it didn't make any sense to me how you had it:)

% create a new figure
h = figure;

% modify some axes properties
xlabel('k * \lambda_D')
ylabel(' omega \omega')
hold on

% save the figure (but nothing has been plotted in this figure?!)
saveas(h,sprintf('omega3%d.png',col));