diff --git a/src/sequential.cpp b/src/sequential.cpp
index c98d255815eb909272aaba638d1b2866dee70f5a..f7933d9fa0bca85c0fd34169cfce7e7b0bcc3ccd 100644
--- a/src/sequential.cpp
+++ b/src/sequential.cpp
@@ -90,7 +90,6 @@ gi::ex multiply_1level( tensor3D_t& T, matrix_int_t& J, int size ) {  // simpler
 				for( b1 = 0 ; b1 < size ; b1++ ){
                     Tn = one_level1_product( &T, &J, A, size, a1, a2, a3, b1 );
                     Tens += Tn;
-                    std::cout << "T" << A << "=" << Tn << ";" << std::endl;
                 }
 			}
 		}
@@ -103,29 +102,42 @@ gi::ex one_level1_product( tensor3D_t* T, matrix_int_t *J, gi::ex A, int size, i
     int b2, b3, c1, c2, c3, d1, d2, d3;
 
     gi::ex Tens = 0;
+    gi::ex T0, T1, T2, T3, T4, T5;
     
     TAB = (*J)[a1][b1]; 
     for( b2 = 0 ; b2 < size ; b2++ ){
         for( b3 = 0 ; b3 < size ; b3++ ){
             TABB =  TAB * A*(*T)[b1][b2][b3];
+            T5 = 0;
             /* Beyond this point, b1 is not used anymore */
             for( c1 = 0 ; c1 < size ; c1++ ){
+                T4 = 0;
                 for( c2 = 0 ; c2 < size ; c2++ ){
                     TABC = TABB * (*J)[a2][c2];
+                    T3 = 0;
                     for( c3 = 0 ; c3 < size ; c3++ ){
-                        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ; 
+                        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
+                        T2 = 0;
                         for( d1 = 0 ; d1 < size ; d1++ ){
                             TABCD = TABCC * (*J)[c1][d1];
+                            T1 = 0;
                             for( d2 = 0 ; d2 < size ; d2++ ){
                                 TABCDD = TABCD * (*J)[b2][d2];
+                                T0 = 0;
                                 for( d3 = 0 ; d3 < size ; d3++ ){
-                                    Tens += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
+                                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                                 }
+                                T1 += T0;
                             }
+                            T2 += T1;
                         }
+                        T3 += T2;
                     }
+                    T4 += T3;
                 }
+                T5 += T4;
             }
+            Tens += T5;
         }
     }
 
@@ -148,21 +160,28 @@ gi::ex one_level2_product( tensor3D_t* T, matrix_int_t *J, gi::ex A, int size, i
     int c3, d1, d2, d3;
 
     gi::ex Tens = 0;
+    gi::ex T0, T1, T2;
     
     TAB = (*J)[a1][b1]; 
     TABB =  TAB * A*(*T)[b1][b2][b3];
     TABC = TABB * (*J)[a2][c2];
     for( c3 = 0 ; c3 < size ; c3++ ){
-        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ; 
+        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
+        T2 = 0;
         for( d1 = 0 ; d1 < size ; d1++ ){
             TABCD = TABCC * (*J)[c1][d1];
+            T1 = 0;
             for( d2 = 0 ; d2 < size ; d2++ ){
                 TABCDD = TABCD * (*J)[b2][d2];
+                T0 = 0;
                 for( d3 = 0 ; d3 < size ; d3++ ){
-                    Tens += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
+                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                 }
+                T1 += T0;
             }
+            T2 += T1;
         }
+        Tens += T2;
     }
     
     return Tens;
@@ -203,19 +222,26 @@ gi::ex two_level1_product( tensor3D_t* T, matrix_int_t *J, gi::ex A, int size, i
     int b2, b3, c1, c2;
 
     gi::ex Tens = 0;
-    
+    gi::ex T0, T1, T2;
+
     TAB = (*J)[a1][b1]; 
     for( b2 = 0 ; b2 < size ; b2++ ){
+        T2 = 0;
         for( b3 = 0 ; b3 < size ; b3++ ){
             TABB =  TAB * A*(*T)[b1][b2][b3];
+            T1 = 0;
             /* Beyond this point, b1 is not used anymore */
             for( c1 = 0 ; c1 < size ; c1++ ){
+                T0 = 0;
                 for( c2 = 0 ; c2 < size ; c2++ ){
                     TABC = TABB * (*J)[a2][c2];
-                    Tens += two_level2_product( T, J, TABC, size, a3, b2, b3, c1, c2 );
+                    T0 += two_level2_product( T, J, TABC, size, a3, b2, b3, c1, c2 );
                 }
+                T1 += T0;
             }
+            T2 += T1;
         }
+        Tens += T2;
     }
     return Tens;
 }
@@ -226,18 +252,25 @@ gi::ex two_level2_product( tensor3D_t* T, matrix_int_t *J, gi::ex TABC, int size
     gi::ex TABCC, TABCD, TABCDD;
     
     gi::ex Tens = 0;
-    
+    gi::ex T0, T1, T2;
+
     for( c3 = 0 ; c3 < size ; c3++ ){
-        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ; 
+        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
+        T2 = 0;
         for( d1 = 0 ; d1 < size ; d1++ ){
             TABCD = TABCC * (*J)[c1][d1];
+            T1 = 0;
             for( d2 = 0 ; d2 < size ; d2++ ){
                 TABCDD = TABCD * (*J)[b2][d2];
+                T0 = 0;
                 for( d3 = 0 ; d3 < size ; d3++ ){
-                    Tens += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
+                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                 }
+                T1 += T0;
             }
+            T2 += T1;
         }
+        Tens += T2;
     } 
     return Tens;
 }