I am having trouble computing the steady state of a model. I can solve it for certain parameter values but not for others (for instance, the steady state can be computed as it is for the mod file posted below, but if I change a_t = 1 in line 59 the steady state can no longer be computed), even though the residuals at the starting values I provide seem to be small and restricted to only two out of 49 equations. I realize that the best way would have been to solve for the steady state by hand, but I am not able to do so, given the highly complex/non linear nature of some of the equations. In the initval block, I have been able to simply the steady state computation to the point where I need to provide guesses for only two endogenous variables, and the rest can be computed recursively.
I have tried other options like homotropy but they dont seem to work. Any suggestions would be welcome. In particular, I would like an opinion on whether the problem is one of computation (so that moving to a more powerful computation platform can be helpful) or something fundamentally related to the model.

The parameter change you are considering seems to be rather large. Are you sure that a steady state exists for the parameter computation? Given the size of the model, I would try to go for a steady state file that uses a solver for the last few equations as is done in the NK_baseline.mod in the Dynare examples folders.

I tried the steady state file approach using fsolve and while its a bit better, I continue to face the problem of computing the steady state for most parameter values.
As mentioned before however, I can boil down the steady state computation to a problem of solving two non linear equations in two unknowns. I wanted an opinion on the following procedure, given that I am not able to compute the exact steady state, even numerically.

1.Instead of using fsolve or something similar, I solve the problem by doing a grid search in 2 dimensions.
2. Step (1) gives me a steady state which is not exact, so lets call it an approximate steady state.
3. Since this is not be an exact steady state, dynare prevents me from running a stochastic simulation based on this. But what if I use a brute force approach here- so that I linearize the model around the approximate steady state and feed the linear model into dynare and run a stochastic simulation to get impulse responses?

I wanted to know if this is sounds even remotely reasonable or are there glaring technical pitfalls here that I am completely overlooking.

I think there are two main issues (perhaps related) that I am still left facing even though the problem is simplified to a 2 equation system.
(a)complex roots, which I frequently encounter
(b) multiple equilibria-its a model of banking with default so a higher interest rate and a higher default rate is also an equilibrium. I think some unreasonable steady states I am getting are due to this issue-the global minimum happens to be at an inefficient point and the algorithm picks it out, whereas the solution that I am looking for is perhaps just a local minimum.