function [q1,q2,dq1,dq2] = DirDin2R_RK4(z1,z2,z3,z4,T1,T2)

global h

%-- Avaliação dos 1's ---------------------------------------------------
 [M1,N1] = InerciasTCG(z1,z2,z3,z4,T1,T2);
 M1=h*M1;  N1=h*N1;  K1=h*z3;  L1=h*z4;
 
%-- Avaliação dos 2's ---------------------------------------------------
 [M2,N2] = InerciasTCG(z1+K1/2,z2+L1/2,z3+M1/2,z4+N1/2,T1,T2);
 M2=h*M2;  N2=h*N2;  K2=h*(z3+M1/2);  L2=h*(z4+N1/2);
 
%-- Avaliação dos 3's ---------------------------------------------------
 [M3,N3] = InerciasTCG(z1+K2/2,z2+L2/2,z3+M2/2,z4+N2/2,T1,T2);
 M3=h*M3;  N3=h*N3;  K3=h*(z3+M2/2);  L3=h*(z4+N2/2);
 
%-- Avaliação dos 4's --------------------------------------------------- 
 [M4,N4] = InerciasTCG(z1+K3,z2+L3,z3+M3,z4+N3,T1,T2);
 M4=h*M4;  N4=h*N4;  K4=h*(z3+M3);  L4=h*(z4+N3);
 
%-- Formulas de Runge-Kutta ---------------------------------------------
 q1 = z1+(K1+2*K2+2*K3+K4)/6.0;
 q2 = z2+(L1+2*L2+2*L3+L4)/6.0;

dq1 = z3+(M1+2*M2+2*M3+M4)/6.0;
dq2 = z4+(N1+2*N2+2*N3+N4)/6.0;
%--------------------------------------------------------------------end
