function dYdt = spring(t,Y)
%  'Y' is a vector of form y = [ r ;  th ; vr ; vth ]

x    = [Y(1);Y(2)];
v    = [Y(3);Y(4)];

global g m l k

a = [((x(1)+l)*(v(2))^2) + (g*cos(x(2))) - (k * x(1)/m); (-2*v(1)*v(2) - g*sin(x(2)))/(l+x(1))];

dYdt = [ v(1) ; v(2) ; a(1) ; a(2)]; 