function xp = cap7_iderobot2gdl(t,x)
global tau1 tau2

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; %parámetros del robot

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;   
   
theta1=m1*lc1*lc1+m2*l1*l1+m2*lc2*lc2+I1+I2;
theta2=l1*m2*lc2; 
theta3=m2*lc2*lc2+I2; 
theta4=g*(lc1*m1+m2*l1); 
theta5=g*m2*lc2; 
theta6=b1; 
theta7=b2; 
theta8=fc1; 
theta9=fc2; 
theta10=fe1; 
theta11=fe2;
   
%matriz de inercia
M=[theta1+2*theta2*cos(q2), theta3+theta2*cos(q2);
   theta3+theta2*cos(q2), theta3];
 
%matriz de Coriolis y fuerzas centrípetas
 C=[ -2*theta2*sin(q2)*qp2, -theta2*sin(q2)*qp2;
     theta2*sin(q2)*qp1, 0];

 
gq11=theta4*sin(q1)+theta5*sin(q1+q2); gq21=theta5*sin(q1+q2);
 gq=[gq11; gq21]; %par gravitacional
 
%par de fricción viscosa, Coulomb y estática
 fr=[theta6*qp1+theta8*sign(qp1)+theta10*(1-abs(sign(qp1)));
     theta7*qp2+theta9*sign(qp2)+theta11*(1-abs(sign(qp2)))]; 
           
    
tau1=(1-exp(-0.8*t))*29.0+ 68*sin(16*t+0.1) + 9*sin(20*t+0.15);
tau2=(1-exp(-1.8*t))*1.2+ 8*sin(26*t+0.08)+2*sin(12*t+0.34);    

tau=[tau1; tau2]; %señal de excitación persistente
qpp = M^(-1)*(tau-C*qp-gq-fr); %aceleración articular

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