Piston problem (main program source for PC and Mac workstation)

%
%   This program is for 3.175 in "An Introduction to Dynamics, 3rd ed."
%   By David J. McGill and Wilton W. King.
%
%   part 1. constant torque (moment) on bar1
%

%
%   necessary input constant 
%

m = 1;                      %  mass of the piston
l1 = 1;                     %  length of the bar1 (left bar)
l2 = 2;                     %  length of the bar2 (right bar)
M = 5;                      %  constant moment at bar1
I1 = 3;                     %  inertia of bar1

t_final = 10;               %  time duration for calculation
                            %  initial time is always 0 

%
%   Intitial conditions for theta and theta_dot
%

theta = 0;                  %  initial displacement of the bar
theta_dot = 0;              %  initial angular velocity of the bar1
theta_ddot = M/I1;          %  initial angular acceleration of the bar1

%
%   Setup integration solver (ODE45)
%

tol = 1e-20;                %  tolerance for ODE45 solver
                            %  smaller it is, more accurate the solution
                            %  will be
x0 = [theta,theta_dot,theta_ddot];        


%
%  Perform ODE45 solver
%

[t,x] = ode45('p3_175s1',0,t_final,x0,tol);
[n,m] = size(x);

Theta = x(:,1);
Phi = asin((l1/l2)*sin(Theta));
X = l1*cos(Theta)+l2*cos(Phi);

%
%  Plot the results
%

figure(1); clf; orient tall;
subplot(2,1,1),plot(t,Theta);
title('Problem 3.175'); 
xlabel('Time (sec.)');    
ylabel('Theta');

subplot(2,1,2),plot(t,Phi);
ylabel('Phi');
xlabel('Time (sec.)');      

figure(2); clf; orient tall;
subplot(2,1,1),plot(t,X);
title('Problem 3.175');
xlabel('Time (sec.)');
ylabel('X (m)');

subplot(2,1,2),plot(Theta,X);  
xlabel('Theta (rad.)');
ylabel('X (m)');

mxx = max(X);       mnx = min(X);
mxt = max(Theta);   mnt = min(Theta);
axis([mnt,mxt,mnx,mxx])