Shear force, bending moment and elastic curve: Matlab program for cantilever beam

Problem 6.31, Hibbeler's "Mechanics of Materials", 8th. edition.

Solution of a cantilever beam, supported on its left end.
___________________________________________________________________

% hibbeler_6_31.m

% Luis Adan Villa Chaparro / Ricardo E. Avila, 24 Feb. 2014
% 

% Problem 6.31 of Hibbeler's "Mechanics of Materials", 8th. edition.

% A cantilever beam of unit length, supported only at its left end,
% has a uniformly distributed unit load w0 on the left half
% of its length; then the load decreases uniformly from w0 to 
% become zero load on the right half of the beam.

clc

clear
disp('Problem 6.31, Hibbeler''s "Mechanics of Materials"')
disp(' ')

L = 1;      % The length of the beam is set to be unit length

w0 = 1;     % The unit load, which may be increased at will.
dx = L/100; % Division of the beam length into segments.
x = [0 : dx : L]';  % Generate values of the x-axis.

n = size(x, 1);     % Size of the x-axis column of data

V = zeros(n, 1);    % Initialize column of shear force data.
M = V;              % Initialize column of bending moment data.
Y = V;              % Initialize column of elastic curve data.

% Processing of the data: shear force, bending moment

% and elastic curve equation, obtained by double integration.
% The two constants of integration are found to be zero, for a
% cantilever beam supported on its left end, at x = 0.

for index = 1 : n

  % First half of the beam length
  V(index) = .75*w0*L - w0*x(index);
  M(index) = (-7*w0*L^2)/24 + .75*w0*L*x(index) - .5*w0*x(index)^2;
  Y(index) = (-7*w0*L^2*(x(index))^2)/48 + ...
             (3*w0*L*(x(index))^3)/24 - (w0*x(index)^4)/24;

  % Second half of the beam length:
  length_2 = x(index) - L/2;  % Define this variable to save time
  if x(index) >= L/2
      V(index) = V(index) + w0*length_2^2/L;
      M(index) = M(index) + w0*length_2^3/(3*L);
      Y(index) = Y(index) + w0*length_2^5/(60*L) ;
  end
end

E = 1; % Elastic constant is unit value

I = 1; % Moment of Inertia is unit value

figure(1)

plot(x, V, 'r', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('V(x), force units')
title('Shear Force Diagram')

figure(2)

plot(x, M, 'b', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('M(x), force times distance')
title('Bending Moment Diagram')

figure(3)

plot(x, Y, 'b', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('y(x), length units')
title('Elastic Curve Diagram')

disp('________________________________ ')

disp('Successful run of Matlab program ')
disp(' ')

% end of Matlab program

Figure(1)

Figure(2)

Figure(3)