dr11_dt = gt - 2*(gp + gt)*r11 - E1*i*r13 + E1c*i*r31 + G*r33;
dr12_dt = (-gp - 2*gt - d1*i + da*i)*r12 - Ef*i*r13 + E1c*i*r32;
dr13_dt = -(E1c*i*r11) - Efc*i*r12 + (-G - gp - 2*gt - d1*i)*r13 + E1c*i*r33;
dr22_dt = gt - 2*gt*r22 - Ef*i*r23 + Efc*i*r32 + G*r33;
dr23_dt = -(E1c*i*r21) - Efc*i*r22 + (-G - 2*gt - da*i)*r23 + Efc*i*r33;
dr33_dt = 2*gp*r11 + E1*i*r13 + Ef*i*r23 - E1c*i*r31 - Efc*i*r32 - 2*(G + gt)*r33;