%Program to solve toy beam bending 
%problem
%
%sigma_x calculated using vector
%operations

%Calculate reactions
A = [1 1; 0 12];
B = [400; 2400];
R = A\B; %Reaction vector

%Calculate sigma_x
M = -600;
ro = 2e-2;
ri = 1e-2*(0.5:0.1:1.5);
I = pi*(ro^4 - ri.^4)/4;
sigma_x = 1e-6*M*ro./I

%Plot sigma_x vs. ri
figure(1); %Create figure #1
clf; %Clear current figure
h=plot(1e2*ri,sigma_x,'-r');
set(gca,'Box','on','LineWidth',2,...
    'FontName','Helvetica',...
	'FontSize',14); 
    %Set axis properties
set(h,'LineWidth',2); 
    %Set linewidth for curve
xlabel('r_i (cm)'); 
ylabel('\sigma_x at O (MPa)');
title('Bending stress');
axis square; %Make axis box square