function res = spring_pendulum_der(t,y)

g = 1; % gravity - not needed in the derivative
k = 100; % spring force
l = 5; % length of spring
d = 1000; % damping

ny = sqrt(y(1)^2 + y(2)^2);
sy = y(1)*y(3)+y(2)*y(4);

res = zeros(4,4);

res(1,3) = 1;
res(2,4) = 1;


res(3,1) = -k*(1-l/ny)-k*l*y(1)^2/ny^3 - d*(sy+y(1)*y(3))/ny^2 + 2*d*y(1)^2*sy/ny^4;
res(3,2) = -k*l*y(1)*y(2)/ny^3 - d*y(4)*y(1)/ny^2 + 2*d*y(1)*y(2)*sy/ny^4;
res(3,3) = -d*y(1)^2/(ny^2);
res(3,4) = -d*y(1)*y(2)/(ny^2);

res(4,1) = -k*l*y(2)*y(1)/ny^3 - d*y(3)*y(2)/ny^2 + 2*d*sy*y(2)*y(1)/ny^4;
res(4,2) = -k*(1-l/ny)-k*l*y(2)^2/ny^3 - d*(sy+y(2)*y(4))/ny^2 + 2*d*sy*y(2)^2/ny^4;
res(4,3) = -d*y(1)*y(2)/(ny)^2;
res(4,4) = -d*y(2)^2/(ny)^2;
