Mohr circle using Matlab

This is a proposed solution to problem 1-33 of the textbook
"Mechanics of Materials", 8th. edition, by R. C. Hibbeler,
using a simple Matlab program. The program solves a system
of two linear equations, for each value of the angle, in order to
calculate the values of normal and shear stress; these values are
functions of the angle 'theta'. The angle 'theta' is formed between
a vector, perpendicular to the area on which the stresses are being
observed or calculated, and a horizontal or longitudinal axis in
the material.

Mohr's circle for values of stress is generated by the solution of the problem,
since the normal stress acts as a sine function, and the shear stress
like a cosine function.

The program can be run with unity input values, for lack of other information:
base = 1
height = 1
applied load = 1

% prob1_33.m

% Ricardo E. Avila / University of Ciudad Juarez / Aug. 28, 2013
% ricardo_avila@hotmail.com
% Chihuahua, Mexico
% Problem 1.33, Hibbeler, 'Mechanics of Materials', 8th. ed.

% 'theta' is the angle in radians of a vector normal to

% the inclined area = base * height, A = b * h;
% the 'theta' angle is measured with respect to the longitudinal
% axis of the loaded material.  

clc

clear

disp(' ')

disp('Problem 1.33, Hibbeler, Mechanics of Materials, 8th. ed.')
disp('________________________________________________________')

% Data input

b = input('Base of material dimension: ');
h = input('Height of material dimension: ');
P = input('Applied load: ');

% Data processing

delta_theta = pi/180; % Generate data every 1 degree of angle
theta = (0: delta_theta : pi)';     % Generate column array: theta
sigma = zeros(size(theta,1), 1);    % sigma: normal stress
tau = sigma;                        % tau: shear stress
A = zeros(2, 2);                    % 2x2 matrix
RHS = zeros(2, 1);                  % Right-hand side of linear system
RHS(1,1) = P;                       % Applied load
bh = b * h;                         % Material: base * height

% Solve linear system for solution of equations

for index = 1: size(theta,1)
  A(1, 1) = bh ;
  A(1, 2) = bh * tan(theta(index));
  A(2, 1) = tan(theta(index));
  A(2, 2) = -1;
  stress = A \ RHS;             % Solve linear system of equations
  sigma(index) = stress(1, 1);  % normal stress
  tau(index) = stress(2, 1);    % shear stress
end

% Post-processing of data

figure(1)
plot(theta, sigma)  % Plot variation of normal stress with angle theta
grid
xlabel('theta, radians')
ylabel('sigma, normal stress')

figure(2)

plot(theta, tau)    % Plot variation of shear stress with angle theta
grid
xlabel('theta, radians')
ylabel('tau, shear stress')

figure(3)

% Generate Mohr's circle with values of (sigma, tau)
line(sigma, tau)
axis equal
grid on
xlabel('sigma, normal stress')
ylabel('tau, shear stress')

disp(' ')

disp('*** Successful program execution ***')
disp('____________________________________')

figure(1)

 figure(2)

 figure(3)