Fix for steady

1.) Fixes cases where model was not truly linear, but Dynare did not catch mistake because fvec and jacob contained cancelling NaN
2.) Use correct error code 22 if NaN appear instead of generic 20

The following mod-file from the forum shows the problem:

```
var a b c i k n q r w y;
predetermined_variables k;
varexo epsilon;

parameters alpha beta delta phi phil rho sigma sigmaa Css Iss Wss Yss SIGMA;
alpha=.33;
beta=0.99;
delta=.025;
phi=.6;
phil =1.1;
rho=.5;
sigma=.4;
sigmaa=.005;
SIGMA=(((1/beta)-1+delta)/alpha)^(1/(alpha-1));
Ass=1;
Bss=0;
Nss=((SIGMA^alpha-delta*SIGMA)*(1/(SIGMA^alpha*(1-alpha)))^(1/sigma))^(-sigma/(phi+sigma));
Kss=SIGMA*Nss;
Qss=beta;
Rss=1/beta-1+delta;
Yss=Nss*(Kss/Nss)^alpha;
Wss=SIGMA^alpha*(1-alpha);
Iss=delta*Kss;
Css=(Nss^phi/Wss)^(-1/sigma);

model(linear);

-sigma*c=-sigma*c(+1)-q; //Euler-Gleichung 1
-sigma*c=phi*n-w; //Euler-Gleichung 2
-sigma*c=-sigma*c(+1)+(1-beta*(1-delta))*r(+1); //Euler-Gleichung 3
k(+1)=(1-delta)*k+delta*i; //Kapitalakkumulation

y=a+(1-alpha)*n+alpha*k; //Produktionsfunktion
w=((1-alpha)*(k/n)^alpha*(a+alpha*k+alpha*n)-phil*n+phil*n(-1)+beta*phil*(n(+1)-n))/Wss;
r=a+(1-alpha)*n-(1-alpha)*k;//Gewinnmaximierungsbedinung 2

y=Css/Yss*c+Iss/Yss*i; //Markträumung 1
b=0; //Markträumung 2
a=rho*a(-1)+epsilon; //Technologischer Fortschritt

end;

steady;

check;

shocks;
var epsilon; 
stderr sigmaa;
end;

stoch_simul (drop=0,periods=0, irf=0);
```

Equation 6 is not linear and division by 0 occurs. Thus, one of the residuals is NaN. But the residual test uses a max operator that ignores NaN and returns ys_init, which is a vector of zeros. This makes the problem hard to spot. The additional test for no NaN should not be too computationally burdensome.