s1=1;
J1=1.5e21;
n01=2.6e18;
tau1=7.2e-4;
f1=1.7e24;
V1=3.8e-29;
sigma1=4.95e-21;

s2=1;
J2=1.5e19/2;
n02=1.1e18;
tau2=.00009;
f2=f1;
V2=6e-28;
sigma2=2.3e-18;

kd1=s1*J1/n01+1/tau1+sigma1*f1;
kr1=s1*J1/n01+1/tau1;
kd2=s2*J2/n02+1/tau2+sigma2*f2;
kr2=s2*J2/n02+1/tau2;
kc1=-s1*J1/n02;
kc2=-s2*J2/n01;

nr1=(s1*J1*kr2+s2*J2*kc1)/(kr1*kr2-kc1*kc2);
nr2=(s2*J2*kr1+s1*J1*kc2)/(kr1*kr2-kc1*kc2);
nd1=(s1*J1*kd2+s2*J2*kc1)/(kd1*kd2-kc1*kc2);
nd2=(s2*J2*kd1+s1*J1*kc2)/(kd1*kd2-kc1*kc2);

Dn1=nr1-nd1;
Dn2=nr2-nd2;

kda=(kd1+kd2)/2;
lambda=sqrt((kd1-kd2)^2/4+kc1*kc2);

Xi1=(kc1*Dn2+(kda-kd1)*Dn1)/lambda;
Xi2=(kc2*Dn1+(kda-kd2)*Dn2)/lambda;

ta=10.^[-9:.1:0];

n1hyp=nd1+exp(-kda*ta).*(Dn1*cosh(lambda*ta)+Xi1*sinh(lambda*ta));
n2hyp=nd2+exp(-kda*ta).*(Dn2*cosh(lambda*ta)+Xi2*sinh(lambda*ta));

loglog(ta,n1hyp,'b-',ta,n2hyp,'g-');