diff -ur src/plugins/fits/non_linear.h src/plugins/fits/non_linear.h
--- src/plugins/fits/non_linear.h	2014-02-13 18:41:44.000000000 +0900
+++ src/plugins/fits/non_linear.h	2016-03-04 00:34:31.000000000 +0900
@@ -177,37 +177,45 @@
             }
             iIterations++;
           } while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
-          gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );
-
-          //
-          // determine the fitted values...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            dXInitial[i] = gsl_vector_get( pSolver->x, i );
-          }
-
-          for( i=0; i<iLength; i++ ) {
-            vectorOutYFitted->value()[i] = function_calculate( pInputX[i], dXInitial );
-            vectorOutYResiduals->value()[i] = pInputY[i] - vectorOutYFitted->value()[i];
+            
+            
+          gsl_matrix *J = gsl_matrix_alloc(pSolver->fdf->n, pSolver->fdf->p);
+          if (J != NULL) {
+              gsl_multifit_fdfsolver_jac(pSolver, J);
+                
+              gsl_multifit_covar( J, 0.0, pMatrixCovariance );
+
+              //
+              // determine the fitted values...
+              //
+              for( i=0; i<NUM_PARAMS; i++ ) {
+                dXInitial[i] = gsl_vector_get( pSolver->x, i );
+              }
+
+              for( i=0; i<iLength; i++ ) {
+                vectorOutYFitted->value()[i] = function_calculate( pInputX[i], dXInitial );
+                vectorOutYResiduals->value()[i] = pInputY[i] - vectorOutYFitted->value()[i];
+              }
+
+              //
+              // fill in the parameter values and covariance matrix...
+              //
+              for( i=0; i<NUM_PARAMS; i++ ) {
+                vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
+                for( j=0; j<NUM_PARAMS; j++ ) {
+                  vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+                }
+              }
+
+              //
+              // determine the value of chi^2/nu
+              //
+              scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
+
+              bReturn = true;
+              
+              gsl_matrix_free( J );
           }
-
-          //
-          // fill in the parameter values and covariance matrix...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
-            for( j=0; j<NUM_PARAMS; j++ ) {
-              vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
-            }
-          }
-
-          //
-          // determine the value of chi^2/nu
-          //
-          scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
-
-          bReturn = true;
-
           gsl_matrix_free( pMatrixCovariance );
         }
         gsl_multifit_fdfsolver_free( pSolver );
diff -ur src/plugins/fits/non_linear_weighted.h src/plugins/fits/non_linear_weighted.h
--- src/plugins/fits/non_linear_weighted.h	2014-02-13 18:41:44.000000000 +0900
+++ src/plugins/fits/non_linear_weighted.h	2016-03-04 00:37:41.000000000 +0900
@@ -193,37 +193,43 @@
           }
           while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
 
-          gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );
-
-          //
-          // determine the fitted values...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            dXInitial[i] = gsl_vector_get( pSolver->x, i );
-          }
-
-          for( i=0; i<iLength; i++ ) {
-            vectorOutYFitted->value()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
-            vectorOutYResiduals->value()[i] = pInputs[YVALUES][i] - vectorOutYFitted->value()[i];
-          }
-
-          //
-          // fill in the parameter values and covariance matrix...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
-            for( j=0; j<NUM_PARAMS; j++ ) {
-              vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
-            }
+          gsl_matrix *J = gsl_matrix_alloc(pSolver->fdf->n, pSolver->fdf->p);
+          if (J != NULL) {
+              gsl_multifit_fdfsolver_jac(pSolver, J);
+                
+              gsl_multifit_covar( J, 0.0, pMatrixCovariance );
+
+              //
+              // determine the fitted values...
+              //
+              for( i=0; i<NUM_PARAMS; i++ ) {
+                dXInitial[i] = gsl_vector_get( pSolver->x, i );
+              }
+
+              for( i=0; i<iLength; i++ ) {
+                vectorOutYFitted->value()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
+                vectorOutYResiduals->value()[i] = pInputs[YVALUES][i] - vectorOutYFitted->value()[i];
+              }
+
+              //
+              // fill in the parameter values and covariance matrix...
+              //
+              for( i=0; i<NUM_PARAMS; i++ ) {
+                vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
+                for( j=0; j<NUM_PARAMS; j++ ) {
+                  vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+                }
+              }
+
+              //
+              // determine the value of chi^2/nu
+              //
+              scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
+
+              bReturn = true;
+              
+              gsl_matrix_free( J );
           }
-
-          //
-          // determine the value of chi^2/nu
-          //
-          scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
-
-          bReturn = true;
-
           gsl_matrix_free( pMatrixCovariance );
         }
         gsl_multifit_fdfsolver_free( pSolver );