How to make a shear force diagram (SFD) and bending moment diagram (BMD) using Matlab?

I am very confused on how to use a written problem and convert it in order to find the SFD and the BMD using Matlab. I tried looking at many Youtube videos, but coudln't figure it out.

This is the problem I used:

=

= -----EQ.1

= (input into EQ.1)

Ans:

To find x:

36.5kN=8(x)

x = 4.56

Shear Force Diagram:

AB:

BC:

Bending Moment Diagram:

AB:

BC:

This is the program I am using. My professor want us to use this program.

function main
  
  %first, we generate the domain and number of points for plotting
  xmin = 0;  %Starting Length of the beam?%
  xmax = 8; %Length of the beam?%
  npt = 4.5625;  
  x = linspace(xmin, xmax, npt);
  A1  = 39.5; 
  C1 = 36.5;
  
  %NOW, ENTER YOUR DESIRED M(x) RELATION HERE:
  %----------------------------------------------------
  V = A1-8(x);
  M = A1(x)-8(x)(x/2);
  %----------------------------------------------------
    
  % the following to the rest will automatically plot your M and V
  figure(1);
  plot(x, M, 'k--');
  xlabel('x');
  ylabel('M(x)');
  figure(2);
  V = diff(M)./diff(x);
  plot(x(2:end), V, 'b-.');
  xlabel('x');
  ylabel('V(x)');
  
  % done
end

%The following function, tells the dumb computer about our structural 
% definition of the Heaviside function 
function out = hevi(x)
  out = zeros(size(x));
  for i=1:length(x)
    if ( x(i) >= 0.)
      out(i) = x(i);
    else
      out(i) = 0.;
    end
  end
end

Answer 1

[P]

clear ALL
clc
disp('Simply Supported Beam');
disp(' ');
L = input('Length of beam in meter = ');
disp(' ');disp('Type 1 for point load, Type 2 for udl')
Type = input('Load case = ');

if Type == 1
disp(' ');
W = input('Load applied in kN = ');
disp(' ');
a = input('Location of Load from left end of the beam in meter = ');
c = L-a;

R1 = W*(L-a)/L;
R2 = W*a/L;

else
disp(' ');
W = input('Uniformly distributed load in kN/m = ');
disp(' ');
b = input('Length of udl in meter = ');
disp(' ');
cg = input('C.G of udl from left end of the beam in meter = ');
a = (cg-b/2);
c = L-a-b;

R1 = W*b*(b+2*c)/(2*L); 
R2 = W*b*(b+2*a)/(2*L); 
end


n = 1000;
delta_x = L/n;
x = (0:delta_x:L)'; 

V = zeros(size(x, 1), 1); 
M = zeros(size(x, 1), 1); 


if Type == 1
for ii = 1:n+1

V(ii) = R1;
M(ii) = R1*x(ii);


if x(ii) >= a
V(ii) = R1-W;
M(ii) = R1*x(ii)-W*(x(ii)-a);
end
end
x1 = a;
Mmax = W*a*(L-a)/L;
else
for ii = 1:n+1

if x(ii) < a
V(ii) = R1;
M(ii) = R1*x(ii);
elseif a <= x(ii) && x(ii)< a+b

V(ii) = R1-W*(x(ii)-a);
M(ii) = R1*x(ii)-W*((x(ii)-a)^2)/2;
elseif x(ii) >= (a+b)

V(ii) = -R2;
M(ii) = R2*(L-x(ii));
end
end
x1 = a+b*(b+2*c)/(2*L);
Mmax = W*b*(b+2*c)*(4*a*L+2*b*c+b^2)/(8*L^2);
end 
subplot(2,1,1);
plot(x, V, 'b-', 'LineWidth', 2);
title('Shear Force Diagram');
xlabel('Position (m)');
ylabel('Shear Force (N)');
grid on;

subplot(2,1,2);
plot(x, M, 'r-', 'LineWidth', 2);
title('Bending Moment Diagram');
xlabel('Position (m)');
ylabel('Bending Moment (N.m)');
grid on;