function [t , X ] = spring_pendulum_int(tspan,x0,l_0,m,k )
%set constants and defaults
if ~exist('l_0','var') , l_0 = 0.5;end
if ~exist('m','var') , m = 1;end
if ~exist('k','var') , k = 5;end
g = 9.81; %acceleration due to gravity m/s^2
%set initial conditions and time span
if ~exist('x0','var') , x0 = [1 ,0 ,pi/4 ,0];end
if ~exist('tspan','var') , tspan = [0 ,20];end
%integrate equations of motion
[t , X ] = ode45( @springpendulum_eq , tspan , x0 );
%integrator for equations of motion
function dx = springpendulum_eq(t , x )
%y = [r, rd, th, thd]
r = x(1);
rd = x(2);
th = x(3);
thd = x(4);
dx = [ rd;
r * thd ^2 + g *cos( th ) - k / m *( r - l_0 );
thd ;
-2* rd * thd / r - g *sin( th )/ r ];
end
end