ode3.m

function dr=ode3(t,r,l,k,m)
  % first order ODE for pendulum attached to 2 springs 
  % l is pendulum length
  % k is spring stiffness of one spring
  % m is mass of pendulum bob
  g=9.81;
  dr=zeros(size(r));
  dr(1)=r(2); % dtheta/dt=dtheta/dt
  dr(3)=r(4); % dy/dt=dy/dt
  if abs(r(1))<=0.2
    dr(2)=-2*k*g/(2*k*l+m*g)*r(1);
    dr(4)=-2*k/m*r(3);
  else
    ddy=1/l*1/sin(r(1))^2*[1 -cos(r(1));-l*cos(r(1)) l]*[-g*sin(r(1));-2*k/m*r(3)];
    dr(2)=ddy(1);
    dr(4)=ddy(2);
end