During Bayesian estimation of a search and matching model, if I have 8 variable steady states, (some parameters in which are to be estimated),and also 8 equations to pin down all of those steady states, then generally it is very hard to have analytical solution for all the steady states using pencil and paper. Is there any code than can help to just keep them unsolved and then run estimation?

Dynare by default tries to solve the resulting simultaneous equation system. The problem is that this brute force approach does not work that well for more complicated models.

Could I ask you further about this? I have attached the mod-file I write for replicating search and match model by Lubik(2009).

In the model part, there are 8 variable steady states (with # sign in front) to be solved and there are 8 equations to pin down those values. Those steady states equations also consist some parameters to be estimated. However, I cannot get analytical solution for these values since it is complicated.

Is there any way to solve this problem, like use fsolve or fzero in 4.4.3 or solve_algo=4 in unstable version? Could I have any example for solving steady states numerically during ESTIMATION?

Why are you using a linearized model when you cannot compute the steady state around which you are linearizing?

That being said, you can always use an external steady state file to call a solver, compute the steady state of the nonlinear model and assign them to respective parameters. An example for calling a solver is the NK_baseline.mod
in the Dynare examples folder.