function [dydt] = hrr(t,y)
global N
dydt(1,1)=y(2);
dydt(2,1)=y(3);
dydt(3,1)=y(4);
h=(1+N)*y(3)+(N*(2+N)/(1+N))*y(1);
hp=(1+N)*y(4)+(N*(2+N)/(1+N))*y(2);
hold0=y(3);
hold8=h^2+4*N*y(2)^2;
hold1=4*N*y(3)+((1+2*N)/(1+N))*h;
hold2=h^2+4*y(2)^2;
hold3=y(2)*(h*hp+4*y(2)*y(3));
hold4=y(3)^2+y(2)*y(4);
hold5=3*h*hp^2+8*hp*y(2)*y(3)+4*h*hold4;
hold6=h/(h^2+4*y(2)^2);
hold7=(h*hp+4*y(2)*y(3))^2;
dydt(4,1)= (N*(2+N)/(1+N)^2)*hold0;
dydt(4,1)= dydt(4,1) + (N/hold8)*((1/(1+N)^2)*hold1*hold2 + ...
   (4*(1-N)/(1+N)^2)*hold3 + ((1-N)/N/(1+N))*hold5 + ...
   ((1-N)*(1-3*N)/N^2/(1+N))*hold6*hold7);
dydt(4,1)=-dydt(4,1);