From c91e1f895bcc7831b44442e5417f3ac5ddea0e5b Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Mon, 23 May 2022 15:24:54 +0200
Subject: [PATCH 1/2] trust_region.m: trap case where linear combination with
 weight 0 on infinite value returns NaN

---
 matlab/trust_region.m | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/matlab/trust_region.m b/matlab/trust_region.m
index 6b070d62b..3e2396252 100644
--- a/matlab/trust_region.m
+++ b/matlab/trust_region.m
@@ -327,5 +327,9 @@ else
     end
     % Form the appropriate  convex combination of the Gauss-Newton direction and the
     % scaled gradient direction.
-    x = alpha*x + (1.0-alpha)*min(sgnorm, delta)*s;
+    if alpha>0
+        x = alpha*x + (1.0-alpha)*min(sgnorm, delta)*s;
+    else %prevent zero weight on Inf evaluating to NaN
+        x = min(sgnorm, delta)*s;
+    end
 end
\ No newline at end of file

From c9fd266cbba7e2e8a27ea409609b257e90d92de8 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Wed, 25 May 2022 14:21:28 +0200
Subject: [PATCH 2/2] solve1.m: trap zero Jacobian case

---
 matlab/solve1.m | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/matlab/solve1.m b/matlab/solve1.m
index fa6a3262f..5c2cd1246 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -121,7 +121,20 @@ for its = 1:maxit
     end
     if rcond_fjac < sqrt(eps)
         fjac2=fjac'*fjac;
-        p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
+        temp=max(sum(abs(fjac2)));
+        if temp>0
+            p=-(fjac2+sqrt(nn*eps)*temp*eye(nn))\(fjac'*fvec);
+        else
+            errorflag = true;
+            errorcode = 5;
+            if nargout<3
+                skipline()
+                dprintf('SOLVE: Iteration %s', num2str(its))
+                disp('Zero Jacobian.')
+                skipline()
+            end
+            return
+        end
     else
         p = -fjac\fvec ;
     end