Why the plot is not coming? i want to plot for every Phi value corresponding M_total value.

%% M-Phi Curve
clc;
close all;
% Input data
fc = 20;            % Grade of concrete in N/mm^2
Fe = 415;           % Grade of steel in N/mm^2
P = 300;            % Axial compressive load in kN
%Phi = 0.00004;      % Initial slope in per meter
N = 6;              % Number of bars
D = 500;            % Depth of the section in mm
B = 250;            % Width of the section in mm
e_c = 50;           % Effective cover in mm
d = 16;             % Diameter of the bars in mm

%%
for Phi = 0 : 0.00004
    
    %% Calculate strain under only axial compressive force
    sigma_x = (P*1000)/(D*B);    % Stress under axial compressive force in N/mm^2
    
    % Quadratic equation solving
    c = sigma_x/(0.45*fc);
    a = 1;
    b= -2;
    delta = (b.^2)-(4*a*c);
    if(delta < 0)      
        disp("Delta < 0 The equation does not have a real root");
    elseif (delta == 0)
        disp('The equation has only one real roots');
        disp(-b./(2*a));   
    else
        disp('The equation has two real roots');
        epsilon_1 = 0.002*((-b+sqrt(delta))./(2*a));
        epsilon_2 = 0.002*((-b-sqrt(delta))./(2*a));
    end
    epsilon = min(epsilon_1,epsilon_2);
    disp(epsilon)
    
    %% Finding depth of neutral axis
    epsilon_max = input('Enter maximum stain value :');
    epsilon_min = epsilon_max-D*Phi*0.001;
    xu = (epsilon_max*D)/(epsilon_max-epsilon_min);     % Depth of neutral axis in mm
    if(xu > D)
        disp("Neutral axis outside the section");
    elseif(xu==D)
        disp("Neutral axis is at the edge of the section");
    elseif(0<=xu)&&(xu<=D)
        disp("Neutral axis inside the section");
    else
        disp("Neutral axis is above the section");
    end
    
    %% Concrete stress, strain and force calculation
    n = input('please enter number of strips :');
    Pc_0 = 0;
    Mc_0 = 0;
    for i = 0:n-1
        for j = i
            p_j = (D/n)*j + (D/(2*n))
            epsilon_ci = (epsilon_max/xu)*(xu-p_j)
            sigma_ci = 0.45*fc*(2*(epsilon_ci/0.002)-(epsilon_ci/0.002)^2)
            P_ci = sigma_ci*B*(D/n)*0.001
            if (1 <= (i+1))&&((i+1) <= n/2)
                M_ci = P_ci*(((D/2)-((D/(2*n))*(i+1)))/(i+1))*0.001 
            elseif ((i+1) >= ((n/2)+1))
                M_ci = P_ci*((((D/2)*(i-2))-(D/(2*n))))*0.001*(-1)^i
            end
        end
        Pc_total = Pc_0 + P_ci;
        Pc_0 = Pc_total;
        Mc_total = Mc_0 + M_ci;
        Mc_0 = Mc_total;
    end
    disp(Pc_total)
    disp(Mc_total)
    
    %% Steel stress, strain and force calculation
    s = input('please enter number of bar layers :');
    Ps_0 = 0;
    Ms_0 = 0;
    e_0 = e_c;
    for j = 0:s-1
        e_j = e_0
        epsilon_sj = (epsilon_max/xu)*(xu-e_j)
        if(epsilon_sj<0.00144)
            sigma_sj = (288.7/0.00144)*epsilon_sj
        elseif(0.00144<=epsilon_sj)&&(epsilon_sj<=0.00163)
            sigma_sj = 288.7 + ((306.7-288.7)/(0.00163-0.00144))*epsilon_sj
        elseif(0.00163<=epsilon_sj)&&(epsilon_sj<=0.00192)
            sigma_sj = 306.7 + ((324.8-306.7)/(0.00192-0.00163))*epsilon_sj
        elseif(0.00192<=epsilon_sj)&&(epsilon_sj<=0.00241)
            sigma_sj = 324.8 + ((342.8-324.8)/(0.00241-0.00192))*epsilon_sj
        elseif(0.00241<=epsilon_sj)&&(epsilon_sj<=0.00276)
            sigma_sj = 342.8 + ((351.8-342.8)/(0.00276-0.00241))*epsilon_sj
        elseif(0.00276<=epsilon_sj)&&(epsilon_sj<=0.00380)
            sigma_sj = 351.8 + ((360.9-351.8)/(0.0038-0.00276))*epsilon_sj
        end
        Ps_j = sigma_sj*(N/2)*(3.14/4)*d^2*0.001
        Ms_j = Ps_j*((D/2)-(e_c))*0.001*(-1)^(j)
        e_0 = e_j + ((D-(2*e_c))/(s-1))
        Ps_total = Ps_0 + Ps_j;
        Ps_0 = Ps_total;
        Ms_total = Ms_0 + Ms_j;
        Ms_0 = Ms_total;
    end
    disp(Ps_total)
    disp(Ms_total)
    
    %% Check
    P_total = Pc_total + Ps_total;
    if((P-2)<=P_total)&&(P_total<= (P+2))
        disp('ok')
    elseif(P_total>P)
        disp('Redo with less epsilon_max value')
    elseif(P_total<P)
        disp('Redo with more epsilon_max value')
    end
    M_total = Mc_total + Ms_total;
end
figure
hold on
plot(Phi,M_total)

[Dyuman Joshi]

Format your code properly so that it is easier to understand.

This line is the source of your problem

for Phi = 0 : 0.00004

By default the increament value is 1, and as the next incremented value (i.e. 1) is bigger than the stated value, thus it returns 0.

Phi = 0 : 0.00004

Phi = 0

If you want get values of phi between the range, you need to define an increment value which is smaller than 0.00004

Phi = 0 : 1e-5 : 0.00004

Phi = 1×5

1.0e-04 * 0 0.1000 0.2000 0.3000 0.4000

Note that when phi is 0, xu is Inf. Thus all calculations after that will be affected by this result.

And since you are trying to plot phi vs M_total, you might want to include the plot command inside the for loop command.

[Rajkumar]

after doing the above correction, still it is not showing the graph

[Dyuman Joshi]

Is this the edited code?

https://in.mathworks.com/matlabcentral/answers/1924980-why-in-this-code-for-every-k-value-m_total-value-is-not-changing-and-the-plot-is-not-showing-any-g?s_tid=srchtitle