diff --git a/src/lyapunov.c b/src/lyapunov.c
index 93d38b3..e4c3df6 100644
--- a/src/lyapunov.c
+++ b/src/lyapunov.c
@@ -128,16 +128,29 @@ int lyapunov(
       // extract eigenvalues (diagonal elements of Du_prod)
       for(i=0; i<MATSIZE; i++){
 	lyapunov[i]=log(cabs(Du_prod[i*MATSIZE+i]))/(time-prevtime);
-	// add to average
-	lyapunov_avg[i]=lyapunov_avg[i]*prevtime/time+lyapunov[i]*(time-prevtime)/time;
+      }
+
+      // sort lyapunov exponents
+      qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
+
+      // average lyapunov
+      for(i=0; i<MATSIZE; i++){
+	// exclude inf
+	if(! isinf(lyapunov[i])){
+	  lyapunov_avg[i]=lyapunov_avg[i]*prevtime/time+lyapunov[i]*(time-prevtime)/time;
+	}
       }
 
       // print
-      fprintf(stderr,"% .15e\n",time);
       for(i=0; i<MATSIZE; i++){
 	printf("% .15e % .15e % .15e\n",time, lyapunov[i], lyapunov_avg[i]);
       }
       printf("\n");
+      fprintf(stderr,"% .15e",time);
+      // print largest and smallest lyapunov exponent
+      if(MATSIZE>0){
+	fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
+      }
 
       // set Du_prod to Q:
       LAPACKE_zungrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2);
@@ -151,6 +164,17 @@ int lyapunov(
   return(0);
 }
 
+// comparison function for qsort
+int compare_double(const void* x, const void* y) {
+  if (*(const double*)x<*(const double*)y) {
+    return(-1);
+  } else if (*(const double*)x>*(const double*)y) {
+    return(+1);
+  } else {
+    return(0);
+  }
+}
+
 // Jacobian of u_n -> u_{n+1}
 int ns_D_step(
   _Complex double* out,
diff --git a/src/lyapunov.h b/src/lyapunov.h
index 9d039a3..fb256c3 100644
--- a/src/lyapunov.h
+++ b/src/lyapunov.h
@@ -22,6 +22,9 @@ limitations under the License.
 // compute Lyapunov exponents
 int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, double nu, double D_epsilon, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double starting_time, unsigned int nthreads);
 
+// comparison function for qsort
+int compare_double(const void* x, const void* y);
+
 // Jacobian of step
 int ns_D_step( _Complex double* out, _Complex double* un, _Complex double* un_next, int K1, int K2, int N1, int N2, double nu, double epsilon, double delta, unsigned int algorithm, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, bool irreversible, bool keep_en_cst, double target_en);