function dy = delta(y)
	global m l g;
	theta_1 = y(1);
	theta_2 = y(2);
	p_1		= y(3);
	p_2		= y(4);
	
	d_theta_1 = (6/(m*l^2)) * (2 * p_1 - 3*cos(theta_1 - theta_2) * p_2)/(16-9*cos(theta_1 - theta_2)^2);
	d_theta_2 = (6/(m*l^2)) * (8 * p_2 - 3*cos(theta_1 - theta_2) * p_1)/(16-9*cos(theta_1 - theta_2)^2);
	d_p_1 = -1/2*m*l^2 * ( d_theta_1 * d_theta_2 * sin(theta_1 - theta_2) + 3 * (g/l) * sin (theta_1) );
	d_p_2 = -1/2*m*l^2 * (-d_theta_1 * d_theta_2 * sin(theta_1 - theta_2) +     (g/l) * sin (theta_2) );
	
	dy = [	d_theta_1;
			d_theta_2;
			d_p_1;
			d_p_2 ];