function sprdemo(m,R,b,k,tinit,tdinit,pinit,pdinit)

% If we were not supplied the proper number of arguments we'll
% use some default values
if nargin<8
    m=0.001;       % Mass of 1 gram
    R=0.05;        % Length of 5cm
    b=0.025;       % Spring distance of 2.5 cm
    k=3;           % Spring constant of 3
    tinit=pi/4;    % Theta value of pi/4
    tdinit=0;      % Theta dot of 0
    pinit=0;       % Phi of 0
    pdinit=10*pi;  % Phi dot of 10*pi (5 rev/s)
end

% Close all open figure windows
close all

% Some constants
g=9.8;
I=m*R^2/3;
h=pdinit*sin(tinit)^2;

% Generate 1000 equally spaced values for theta
% ranging from 0 to pi/2
theta=linspace(0,pi/2,1000);

% Calculate the three different components of the
% potential energy equation
cent=2*I*h^2./sin(theta);
grav=2*m*g*R*cos(theta);
spring=4*k*b^2*sin(theta).^2;

% Plot the three components of potential energy and
% the total potential energy
plot(theta, cent, ...
    theta, grav, ...
    theta, spring, ...
    theta, cent + grav + spring);

% Provide a legend so we can tell which graph is which
legend('Centrifugal', 'Gravitational', 'Spring', 'Total');

% The centrifugal term is coming down from +inf and can
% really mess with the range of the axis, in general all
% of the interesting behavior will occur between 0 and
% just over the maximum height of the spring graph
top=max(spring);
axis([0, pi/2, 0, top+1/4*top]);

%Create a new figure window
figure

% Use the ode45 numerical routine to calculate values
% for theta, thetadot, phi, and phidot on the time
% interval 0 to 1 seconds.
[t,x]=ode45(@F, ...
      [0,1], ...
      [pinit,pdinit,tinit,tdinit], ...
      [], ...
      m, R, k, b, g, I);
  
% Plot theta and thetadot versus time
plot(t,x(:,3),t,x(:,4));

% Provide a legend so we can tell which graph is which
legend('Theta', 'Thetadot');

% Write the data returned from ode45 to a text file
fid=fopen('odeoutput.txt', 'w+');
for c = 1:size(t)
    fprintf(fid, '%5f %5f %5f\n', [t(c); x(c,3); x(c,1)]);
end
fclose(fid);


% This function calculates values for x1', x2', x3',
% and x4' using the values for x1, x2, x3, and x4
% that are passed to it in the vector x by the ode45
% routine
function xprime=F(t, x, m, R, k, b, g, I)

% These are just the coefficents of the spring and
% gravitational terms, I've broken them out here in
% an attempt to make the lines below easier to read
q=m*g*R/(2*I);
w=2*k*b^2/I;

% Calculate the derivatives
xprime=[x(2); ...
        ...
        -2*x(4)*x(2)*cos(x(3))/sin(x(3)); ...
        ...
        x(4); ...
        ...
        q*sin(x(3))-(w-x(2)^2)*sin(x(3))*cos(x(3))];