dynare++: add steady state to first part of simulation when

interrupted by non finite result

dynare++/kord/decision_rule.hweb

@@ -300,7 +300,7 @@ TwoDMatrix* simulate(emethod em, int np, const Vector& ystart,
 	@<initialize vectors and subvectors for simulation@>;
 	@<perform the first step of simulation@>;
 	@<perform all other steps of simulations@>;
-	@<add the steady state to all columns of |res|@>;
+	@<add the steady state to columns of |res|@>;
 	return res;
 }
 
@@ -327,10 +327,11 @@ evaluate the polynomial.
 	eval(em, out, dyu);
 
 @ Also clear. If the result at some period is not finite, we pad the
-rest of the matrix with zeros and return immediatelly.
+rest of the matrix with zeros.
 
 @<perform all other steps of simulations@>=
-	for (int i = 1; i < np; i++) {
+        int i=1;
+	while (i < np) {
 		ConstVector ym(*res, i-1);
 		ConstVector dym(ym, ypart.nstat, ypart.nys());
 		dy = dym;
@@ -342,14 +343,16 @@ rest of the matrix with zeros and return immediatelly.
 				TwoDMatrix rest(*res, i+1, np-i-1);
 				rest.zeros();
 			}
-			return res;
+			break;
 		}
+		i++;
 	}
 
-@ Even clearer.
-@<add the steady state to all columns of |res|@>=
-	for (int i = 0; i < res->ncols(); i++) {
-		Vector col(*res, i);
+@ Even clearer. We add the steady state to the numbers computed above and leave the padded columns to zero.
+
+@<add the steady state to columns of |res|@>=
+	for (int j = 0; j < i; j++) {
+		Vector col(*res, j);
 		col.add(1.0, ysteady);
 	}