building a graph( small project)

hello,

I was wondering if someone could help me with a problem I have been given. It is about deflection of beams in a course in college called numerical analysis.

Given: 𝑞0 = 1𝑁𝑚𝑚, E = 140MPa, L = 1000mm, v = -0.7mm.

1.

Draw a graph of x in front of I for 100∙ 10^6 mm^4 ≤ I ≤ 450∙10^6 mm^4. What conclusion can .be physically deduced from the graph?

Now, here is the code for x, I found it numerically (function beam1).I used the newton raphson method.

function X_coordinate = beam1(x0,q0,L,I,E)

epsilon = 1e-4;
flag = 0;
counter = 0;
x = x0;
while flag == 0
    %Check for infinite loop
    counter = counter + 1
    if counter > 10000;
        error('Too many iterations');
    end
    v = -(q0*L)/(3*pi^4*E*I)*(48*L^3*cos((pi*x)/(2*L))-48*L^3+3*pi^3*L*x^2-pi^3*x^3)
    v_d = -(q0*L)/(3*pi^4*E*I)*(-48*L^3*(pi/(2*L))*sin((pi*x)/(2*L))+6*pi^3*L*x-3*pi^3*x^2)
     x = x-v/v_d %solution 
    
    if abs(v/v_d) < epsilon || abs(v) < epsilon 
        flag = 1;
    end
end
X_coordinate = x;
if X_coordinate <=0
    error('No solution found,Please enter valid value')
else
    X_coordinate = x
end

Here is the code for section 1. I need to use the solution of the beam1 function for section 1 numerically, of course, and I have no idea how to proceed from here.

function x_I = MomentI(x0,q0,L,v,E)

epsilon = 1e-4;
flag = 0;
counter = 0;
x = x0
I_x = 1:1000
x_x = 1:1000

while flag == 0
    %Check for infinite loop
    counter = counter + 1
    if counter > 10000;
        error('Too many iterations');
    end
    %if I_x >= 100*10^6 & I_x<=450*10^6
         for x=1:1000
            I_x(x) = -(q0*L)/(3*pi.^4*E*v)*(48*L^3*cos((pi*x)/(2*L))-48*L.^3+3*pi.^3*L*x.^2-pi^3*x.^3)
            I_d(x) = -(q0*L)/(3*pi.^4*E*v)*(-48*L.^3*(pi/(2*L))*sin((pi*x)./(2*L))+6*pi^3*L*x-3*pi.^3*x.^2)
            x_x(x) = x_x-I_x/I_d %solution 
         end
    if abs(I/I_d) < epsilon || abs(I) < epsilon 
        flag = 1;
    end
    %end
end
x_I = x;
if x_I <=0
    error('No solution found,Please enter valid value')
else
    x_I = x
end
x = 1:1000;
hold on
title ('I as a function of x')
xlabel('x[mm]')
ylabel('I[mm]')
x_I = plot(x,I_x)
hold off
end

Answer 1

[Netanel Malihi]

anyone can help?