diff --git a/unsupported/Eigen/src/NonLinear/fdjac1.h b/unsupported/Eigen/src/NonLinear/fdjac1.h
index b9ec75bd4..e697e7ab6 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac1.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac1.h
@@ -22,68 +22,47 @@ int ei_fdjac1(
 
     eps = ei_sqrt((std::max(epsfcn,epsmch)));
     msum = ml + mu + 1;
-    if (msum < n) {
-        goto L40;
-    }
-
-    /*        computation of dense approximate jacobian. */
-
-    for (j = 0; j < n; ++j) {
-        temp = x[j];
-        h = eps * ei_abs(temp);
-        if (h == 0.)
-            h = eps;
-        x[j] = temp + h;
-        iflag = Functor::f(x, wa1);
-        if (iflag < 0)
-            goto L30;
-        x[j] = temp;
-        for (i = 0; i < n; ++i) {
-            fjac(i,j) = (wa1[i] - fvec[i]) / h;
-            /* L10: */
+    if (msum >= n) {
+        /* computation of dense approximate jacobian. */
+        for (j = 0; j < n; ++j) {
+            temp = x[j];
+            h = eps * ei_abs(temp);
+            if (h == 0.)
+                h = eps;
+            x[j] = temp + h;
+            iflag = Functor::f(x, wa1);
+            if (iflag < 0)
+                return iflag;
+            x[j] = temp;
+            fjac.col(j) = (wa1-fvec)/h;
         }
-        /* L20: */
-    }
-L30:
-    /* goto L110; */
-    return iflag;
-L40:
 
-    /*        computation of banded approximate jacobian. */
-
-    for (k = 0; k < msum; ++k) {
-        for (j = k; msum< 0 ? j > n: j < n; j += msum) {
-            wa2[j] = x[j];
-            h = eps * ei_abs(wa2[j]);
-            if (h == 0.) h = eps;
-            x[j] = wa2[j] + h;
-            /* L60: */
-        }
-        iflag = Functor::f(x, wa1);
-        if (iflag < 0) {
-            /* goto L100; */
-            return iflag;
-        }
-        for (j = k; msum< 0 ? j > n: j < n; j += msum) {
-            x[j] = wa2[j];
-            h = eps * ei_abs(wa2[j]);
-            if (h == 0.) h = eps;
-            for (i = 0; i < n; ++i) {
-                fjac(i,j) = 0.;
-                if (i >= j - mu && i <= j + ml) {
-                    fjac(i,j) = (wa1[i] - fvec[i]) / h;
-                }
-                /* L70: */
+    }else {
+        /* computation of banded approximate jacobian. */
+        for (k = 0; k < msum; ++k) {
+            for (j = k; msum< 0 ? j > n: j < n; j += msum) {
+                wa2[j] = x[j];
+                h = eps * ei_abs(wa2[j]);
+                if (h == 0.) h = eps;
+                x[j] = wa2[j] + h;
+            }
+            iflag = Functor::f(x, wa1);
+            if (iflag < 0) {
+                return iflag;
+            }
+            for (j = k; msum< 0 ? j > n: j < n; j += msum) {
+                x[j] = wa2[j];
+                h = eps * ei_abs(wa2[j]);
+                if (h == 0.) h = eps;
+                for (i = 0; i < n; ++i) {
+                    fjac(i,j) = 0.;
+                    if (i >= j - mu && i <= j + ml) {
+                        fjac(i,j) = (wa1[i] - fvec[i]) / h;
+                    }
+                }
             }
-            /* L80: */
         }
-        /* L90: */
     }
-    /* L100: */
-    /* L110: */
     return iflag;
-
-    /*     last card of subroutine fdjac1. */
-
 } /* fdjac1_ */
 
diff --git a/unsupported/Eigen/src/NonLinear/fdjac2.h b/unsupported/Eigen/src/NonLinear/fdjac2.h
index fd6cf1d8b..f2e627d4f 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac2.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac2.h
@@ -8,18 +8,15 @@ int ei_fdjac2(
         Matrix< Scalar, Dynamic, 1 >  &wa)
 {
     /* Local variables */
-    Scalar h;
-    int i, j;
-    Scalar eps, temp;
+    Scalar h, temp;
     int iflag;
 
     /* Function Body */
     const Scalar epsmch = epsilon<Scalar>();
     const int n = x.size();
-    const int m = fvec.size();
+    const Scalar eps = ei_sqrt((std::max(epsfcn,epsmch)));
 
-    eps = ei_sqrt((std::max(epsfcn,epsmch)));
-    for (j = 0; j < n; ++j) {
+    for (int j = 0; j < n; ++j) {
         temp = x[j];
         h = eps * ei_abs(temp);
         if (h == 0.) {
@@ -27,21 +24,11 @@ int ei_fdjac2(
         }
         x[j] = temp + h;
         iflag = Functor::f(x, wa);
-        if (iflag < 0) {
-            /* goto L30; */
+        if (iflag < 0)
             return iflag;
-        }
         x[j] = temp;
-        for (i = 0; i < m; ++i) {
-            fjac(i,j) = (wa[i] - fvec[i]) / h;
-            /* L10: */
-        }
-        /* L20: */
+        fjac.col(j) = (wa-fvec)/h;
     }
-    /* L30: */
     return iflag;
-
-    /*     last card of subroutine fdjac2. */
-
-} /* fdjac2_ */
+}