function yp = F(t,y)

T=1100;

D0=9e+14;
C0=1.5;
K=0.25;
Ef=0.5;
Ed=3;
gam=0.27;
omega=8.989;
k=8.617e-5;
D=D0*exp(-Ed/k/T);
rho=C0*exp(-Ef/k/T);

yp=zeros(2,1);

r = -D*K*rho*(exp(gam*omega/k/T/y(1))-exp(-gam*omega/k/T/y(2))) ...
    /(D*(y(1)+y(2))+K*y(1)*y(2)*log(abs(y(2)/y(1))));

yp(1)=r*y(2);
yp(2)=r*y(1);

