close all
clear all

global m;
global k;
global b;


k = 195000000*(pi*.005^2)  %N/m
b = 2*sqrt(k*m)  % For Critical Damping
%b = b/2;  % Use damping of 1/2 critical value to introduce small transients.


%Solver Options
rtol = 1e-8;
atol = 1e-7*ones(16,1);
options = odeset('RelTol',rtol,'AbsTol',atol);

%Initial Conditions

[T,Y] = ode45(@rhs,[0 30],[0 0 0 0],options);

figure
plot(T,Y(:,1),'-',T,Y(:,2),'-.')
figure
plot(T,Y(:,3),'-',T,Y(:,4),'-.')
F21 = k*(Y(:,3)-Y(:,4));
F12 = -k*(Y(:,3)-Y(:,4));
figure
plot(T,F21,'-')