
function f=dubbel(t,u)

 g=1; m=1;
 r=1; R=2; l=2; 

u1=u(1);
u2=u(2);
u3=u(3);
u4=u(4);

alpha=5*(R-r)-3/2*(R-r)*cos(u1-u3).^2;
beta= l/3-l/10*cos(u1-u3).^2;

f2=1/alpha*(3/2*g*cos(u1-u3).*sin(u3)-l*sin(u1-u3).*u4.^2- ... 
   4*g*sin(u1)-3/2*(R-r)*cos(u1-u3).*sin(u1-u3).*u2.^2);
f4=1/beta*(1/2*(R-r)*sin(u1-u3).*u2.^2-1/2*g*sin(u3)+ ... 
   l/10*sin(u1-u3).*cos(u1-u3).*u4.^2+2/5*g*cos(u1-u3).*sin(u1));

f=[u2;f2;u4;f4];
