function dstate = triple_ode(t,state)
du00 = state(1); du01 = state(2); du02 = state(3);
du10 = state(4); du11 = state(5); du12 = state(6);

du20 = (-2*du10^2*sin(du00 - du01)*cos(du00 - du02)*cos(du01 - du02) - 2*du10^2*sin(du00 - du02)*cos(du00 - du01)*cos(du01 - du02) + 2*du10^2*sin(2*du00 - 2*du01) + du10^2*sin(2*du00 - 2*du02) - du11^2*(-sin(du00 - 3*du01 + 2*du02) + sin(du00 + du01 - 2*du02))/2 - 2*du11^2*sin(du00 - du01)*cos(du01 - du02)^2 + 4*du11^2*sin(du00 - du01) + 2*du11^2*sin(du01 - du02)*cos(du00 - du02) + du12^2*(-sin(du00 - 2*du01 + du02) + sin(du00 + 2*du01 - 3*du02))/4 - du12^2*sin(du00 - du02)*cos(du01 - du02)^2 + 2*du12^2*sin(du00 - du02) - 2*du12^2*sin(du01 - du02)*cos(du00 - du01) - 30*sin(du00)*cos(du01 - du02)^2 + 60*sin(du00) - 40*sin(du01)*cos(du00 - du01) + 20*sin(du01)*cos(du00 - du02)*cos(du01 - du02) + 20*sin(du02)*cos(du00 - du01)*cos(du01 - du02) - 20*sin(du02)*cos(du00 - du02))/(cos(2*du00 - 2*du01) + cos(2*du01 - 2*du02)/2 - 5/2);

du21 = (-2*du10^2*(-sin(-3*du00 + du01 + 2*du02) + sin(du00 + du01 - 2*du02)) + 8*du10^2*sin(du00 - du01)*cos(du00 - du02)^2 - 24*du10^2*sin(du00 - du01) + 12*du10^2*sin(du00 - du02)*cos(du01 - du02) + 8*du11^2*sin(du00 - du01)*cos(du00 - du02)*cos(du01 - du02) - 8*du11^2*sin(2*du00 - 2*du01) - 8*du11^2*sin(du01 - du02)*cos(du00 - du01)*cos(du00 - du02) + 6*du11^2*sin(2*du01 - 2*du02) + du12^2*(-sin(-2*du00 + du01 + du02) + sin(2*du00 + du01 - 3*du02)) - 8*du12^2*sin(du00 - du02)*cos(du00 - du01) - 4*du12^2*sin(du01 - du02)*cos(du00 - du02)^2 + 12*du12^2*sin(du01 - du02) - 240*sin(du00)*cos(du00 - du01) + 120*sin(du00)*cos(du00 - du02)*cos(du01 - du02) - 80*sin(du01)*cos(du00 - du02)^2 + 240*sin(du01) + 80*sin(du02)*cos(du00 - du01)*cos(du00 - du02) - 120*sin(du02)*cos(du01 - du02))/(2*(2*cos(2*du00 - 2*du01) + cos(2*du01 - 2*du02) - 5));

du22 = (du10^2*(sin(-3*du00 + 2*du01 + du02) - sin(du00 - 2*du01 + du02)) + 6*du10^2*sin(du00 - du01)*cos(du01 - du02) + 4*du10^2*sin(du00 - du02)*cos(du00 - du01)^2 - 6*du10^2*sin(du00 - du02) + du11^2*(-sin(-2*du00 + du01 + du02) + sin(2*du00 - 3*du01 + du02)) - 4*du11^2*sin(du00 - du01)*cos(du00 - du02) + 4*du11^2*sin(du01 - du02)*cos(du00 - du01)^2 - 6*du11^2*sin(du01 - du02) + 2*du12^2*sin(du00 - du02)*cos(du00 - du01)*cos(du01 - du02) - du12^2*sin(2*du00 - 2*du02) + 2*du12^2*sin(du01 - du02)*cos(du00 - du01)*cos(du00 - du02) - 3*du12^2*sin(2*du01 - 2*du02)/2 + 60*sin(du00)*cos(du00 - du01)*cos(du01 - du02) - 60*sin(du00)*cos(du00 - du02) + 40*sin(du01)*cos(du00 - du01)*cos(du00 - du02) - 60*sin(du01)*cos(du01 - du02) - 40*sin(du02)*cos(du00 - du01)^2 + 60*sin(du02))/(cos(2*du00 - 2*du01) + cos(2*du01 - 2*du02)/2 - 5/2);

dstate = [ du10; du11; du12; du20; du21; du22 ];
