Urgent: How can i write a code to get graph e Vs P

%CE 405 - CIVIL ENGINEERING PROJECT 2015
%---------------------------------------

clc
clear all
close all

fprintf('        CE 405 - CIVIL ENGINEERING PROJECTS\n')
fprintf('DEVELOPMENT OF LOW COST PRE-STRESSED PRECAST SLAB SYSYTEM\n')
fprintf('\n')
fprintf('\n')
%INPUT PARAMETERS
%----------------

%b1 = Bottom width of the beam section
%b2 = Top width of the beam section
%h1 = Height of the flange
%h2 = Height of the web
%l = length to the neutral axis from bottom of the beam
%e = Eccentricity
%fcu = Concrete grade

%CALCULATION OF THE NEUTRAL AXIS
%-------------------------------
prompt1 = 'What is the Bottom width of the beam section in mm? ';
b1 = input(prompt1)* 10^(-3);
prompt2 = 'What is the Top width of the beam section in mm? ';
b2 = input(prompt2)* 10^(-3);
prompt3 = 'What is the height of the flange in mm? ';
h1 = input(prompt3)* 10^(-3);
prompt4 = 'What is the height of the web in mm? ';
h2 = input(prompt4)* 10^(-3);
prompt8 = 'What is the span of the beam in m? ';
L = input(prompt8);
%prompt5 = 'What is the applied eccentricity of the applied force with respect to neutral axis in mm? ';
%e = input(prompt5)* 10^(-3);
prompt6 = 'What is the unit density of the concrete in kN/m^3? ';
w1 = input(prompt6)*1000;
prompt9 = 'What is the grade of the concrete? ';
fcu = input(prompt9);
prompt7 = 'What is the applied live load to unit weight of the beam in kN/m? ';
w2 = input(prompt7)*1000;
fprintf('\n')
fprintf('\n')

%AREA OF THE SECTION
%-------------------
A = b1*h1+b2*h2;

%CALCULATION OF THE SELF WEIGHT PER UNIT LENGTH OF THE BEAM
%------------------------------------------
W1 = A*w1;

%CALCULATION OF THE BENDING MOVEMENT
%-----------------------------------
Mmax = ((W1+w2)*L^2)/8;
Mmin = (W1*L^2)/8;

%CALCULATIONS OF THE NEUTRAL AXIS
%--------------------------------
Yb = ((h1/2*(h1*b1))+((h2*b2)*(h1+h2/2)))/((h1*b1)+(h2*b2));
Yt = h1+h2-Yb;

% CALCULATION OF SECOND MOMENT AREA
%----------------------------------
I=(b1*(h1^3)/12)+((b1*h1)*(Yt-h1/2)^2)+(b2*(h2^3)/12)+((b2*h2)*(h1+(h2/2)-Yt)^2); %second moment of area

%SECTION MODULUS
%---------------
Zb = I/Yb;
Zt = I/Yt;

PROPERTY= {'Cross section area of the beam'; 'Self weight per unit length of the beam';'Maximum bending moment at middle';'Minimum bending moment at middle';'Distance between neutral axis and bottom of the beam';'Distance between neutral axis and top of the beam';'Second moment of area of the section';'Section modulus of beam section with respect to Yb';'Section modulus of beam section with respect to Yt'};
VALUES = [A; W1; Mmax; Mmin; Yb; Yt; I; Zb; Zt];
data = dataset(PROPERTY, VALUES);
disp(data);

%CALCULATION OF MINIMUM PRE-STRESSED FORCE
%-----------------------------------------

fmax = (fcu*1000000)/3;
fmin = 0;
fmmax = 16*10000000;
fmmin = -1*10000000;

for e=0:0.005:0.1
    a=e;
    P=(fmax-(Mmax/Zt))/(0.8*((1/A)-(a/Zt)));

end

figure
plot(e,P)

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

Answer 1

[Walter Roberson]

e_vals = 0:0.005:0.1;
P = zeros(1,length(e_vals));
for eidx = 1 : length(e_vals)
    a = e_vals(eidx);
    P(eidx)=(fmax-(Mmax/Zt))/(0.8*((1/A)-(a/Zt)));
end

figure
plot(e_vals,P)