Why in this code for every k value, M_total value is not changing? And the plot is not showing any graph?

%% M-Phi Curve

clc;

clear all;

close all;

% Input data

fck = 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

k = 0.00001:0.00001:0.00004;

M_total = zeros(1,length(k));

%% Cross sectional area of the section

A = B*D; % Rectangular section Area in mm^2

%% Calculate strain under only axial compressive force

sigma_x = (P*1000)/(A); % Stress under axial compressive force in N/mm^2

% Quadratic equation solving

c = sigma_x/(0.45*fck);

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 :');

for Phi = k

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 concrete 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*fck*(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 steel bar layers :');

Ps_0 = 0;

Ms_0 = 0;

e_0 = e_c;

fy = input('please enter the grade of steel :');

if (fy == 415)

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.00380-0.00276))*epsilon_sj;

elseif(epsilon_sj>0.00380)

sigma_sj = 360.9;

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

elseif (fy == 500)

for j = 0:s-1

e_j = e_0;

epsilon_sj = (epsilon_max/xu)*(xu-e_j);

if(epsilon_sj<0.00174)

sigma_sj = (347.8/0.00174)*epsilon_sj;

elseif(0.00174<=epsilon_sj)&&(epsilon_sj<=0.00195)

sigma_sj = 347.8 + ((369.6-347.8)/(0.00195-0.00174))*epsilon_sj;

elseif(0.00195<=epsilon_sj)&&(epsilon_sj<=0.00226)

sigma_sj = 369.6 + ((391.3-369.6)/(0.00226-0.00195))*epsilon_sj;

elseif(0.00226<=epsilon_sj)&&(epsilon_sj<=0.00277)

sigma_sj = 391.3 + ((413.0-391.3)/(0.00277-0.00226))*epsilon_sj;

elseif(0.00277<=epsilon_sj)&&(epsilon_sj<=0.00312)

sigma_sj = 413.0 + ((423.9-413.0)/(0.00312-0.00277))*epsilon_sj;

elseif(0.00312<=epsilon_sj)&&(epsilon_sj<=0.00417)

sigma_sj = 423.9 + ((434.8-423.9)/(0.00417-0.00312))*epsilon_sj;

elseif(epsilon_sj>0.00417)

sigma_sj = 434.8;

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

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;

disp(M_total)

disp(P_total)

end

%% Plot

figure

hold on

plot(Phi,M_total)

[Edoardo_a]

I think the problem is that the variable MC_Total is updated everytime inside its for loop.

You should store each time you generate a new MC_Total value somewhere outside the for loop in an array and then plot it against Phi. The plot doesn't show anything because you are plotting a point with coordinate X(Phi) and Y(MC_Total).

If you want to update the figure everytime the for loop in Phi is running, you should put the end of the loop at the very end and the open the figure before you start the loop.