how can I solve equation of motion a pair spur gear with ODE

function xdot=f(t,x)
a=20*(pi/180); % pressure degree
m=0.003; % modul
z1=40; % number of teeth of gear 1
z2=40; % number of teeth of gear 2
Rb1=(m*z1*cos(a))/2; % base circle of gear 1
Rb2=(m*z2*cos(a))/2; % base circle of gear 2
Rv1=m*((z1/2)+1); % tip radius of gear 1
Rv2=m*((z2/2)+1); % tip radius of gear 2
T1T2=(Rb1+Rb2)*tan(a);
T2P1=sqrt(Rv2^2-Rb2^2);
T1P1=T1T2-T2P1;
T1P2=sqrt(Rv1^2-Rb1^2);
T1P=Rb1*tan(a);
T2P=Rb2*tan(a);
T2P2=T1T2-T1P2; 
O1O2=m*(z1+z2)/2;

gama1=a-atan(T1P1/Rb1);
gama2=atan(T2P1/Rb2)-a;
betha1=atan(T1P2/Rb1)-a;
betha2=a-atan(T2P2/Rb2);
n1=1:z1;
n2=1:z2;
theta1=x(1)-n1*gama1-(n1-1)*betha1;
theta2=x(3)+n2*gama2+(n2-1)*betha2;

Rc1=Rb1./cos(a-gama1+theta1);
O1P1=Rb1./cos(a-gama1);
P1Pc=sqrt(O1P1^2+Rc1.^2-2*O1P1*Rc1.*cos(theta1));
PcP1O1=(acos((P1Pc.^2+O1P1^2-Rc1.^2)./(2*P1Pc.*O1P1)));
PcP1O2=pi-gama1-gama2-PcP1O1;
Rc2=sqrt(Rv2^2+P1Pc.^2-2*Rv2*P1Pc.*cos(PcP1O2));
% theta2=-asin(Rc1.*sin(theta1)./Rc2);

v = 0.3; %Poisson’s ratio
E = 200*10^9; %Young’s modulus E
b1 = 0.005;
b2 = 0.005;
hh = 0.001;
tt = 0.01; 
c1 = (E*(1-v))/(b1*(1+v)*(1-2*v));
c2 = (E*(1-v))/(b2*(1+v)*(1-2*v));
c3 = (c1*c2)/(c1+c2);
k = c3*tt*hh; %stiffness coefficient
c = 1000; %damping coefficient
m1 = pi*Rv1^2*tt*7850; % kg
m2 = pi*Rv2^2*tt*7850;
I1 = (1/2)*m1*Rv1^2; 
I2 = (1/2)*m2*Rv2^2;
T1 = 0; %torque
T2 = 0;

xdot = zeros(4,1);
xdot(1) = x(2);
xdot(2) = (T1-(c*(Rc1.*x(2)-Rc2.*x(4))+k*(Rc1.*x(1)-Rc2.*x(3))).*Rc1)/I1;
xdot(3) = x(4);
xdot(4) = (-T2+(c*(Rc1.*x(2)-Rc2.*x(4))+k*(Rc1.*x(1)-Rc2.*x(3))).*Rc2)/I2;
end

i have this fuction, but Rc1 and Rc2 are variable and they've shuold getting change while x does. but i dont know how can i manage this?!

[Torsten]

You will have to define 4*n1*n2 equations instead of 4.