function xp =robot2gdl(t,x)
q1=x(1); q2=x(2); q = [q1; q2]; %vector de posición articular
qp1=x(3); qp2=x(4); qp = [qp1; qp2]; %vector de velocidad articular

m1=23.902; l1=0.45; lc1=0.091; I1=1.266; 
b1=2.288; fc1=7.17; fe1=8.8;
m2=3.880; l2=0.45; 
lc2=0.048; I2=0.093; 
b2=0.175; fc2=1.734; fe2=1.87; g=9.81;   
   
 m11=m1*lc1*lc1+m2*l1*l1+m2*lc2*lc2+2*m2*l1*lc2*cos(q2)+I1+I2; 
 m12=m2*lc2*lc2+m2*l1*lc2*cos(q2)+I2;
 m21=m12;
 m22=m2*lc2*lc2+I2;
 
 M=[m11, m12; m21, m22];
 
 c11=-2*m2*l1*lc2*sin(q2)*qp2;
 c12=-m2*l1*lc2*sin(q2)*qp2;
 c21=m2*l1*lc2*sin(q2)*qp1;
 c22=0;
 
 C=[ c11, c12; c21, c22];
 
 gq11=(lc1*m1++m2*l1)*sin(q1)+m2*lc2*sin(q1+q2);
 gq21=m2*lc2*sin(q1+q2);
 gq=g*[gq11; gq21];
 
 fr=[b1*qp1+fc1*sign(qp1)+fe1*(1-abs(sign(qp1)));
     b2*qp2+fc2*sign(qp2)+fe2*(1-abs(sign(qp2)))]; 
           
    
tau=[(1-exp(-0.8*t))*29.0+ 68*sin(16*t+0.1) + 9*sin(20*t+0.15);
      (1-exp(-1.8*t))*1.2+ 8*sin(26*t+0.08)+2*sin(12*t+0.34) ];    

qpp = M^(-1)*(tau-C*qp-gq-fr);

% Vector de salida
xp = [qp1; qp2; qpp(1); qpp(2)];
end