If my model has 2 parameters , say alpha and beta jointly determine the in/determinacy regions, and I would like to see region(alpha=0,1,2…100 ; beta=0,1,2…5), a slow method is to run stoch_simul or simul and each time calibrate a pair of values of alpha and beta to see if the model is indeterminate or not. In such case I need run 500 times of stoch_simul or simul.

Is there any easier way to do this ? Could you please give me some hints?

You are running into a problem I cannot tell you how many times I have run into. The best advice I can give you is the following:

Use a non-linear system solver (can be done within matlab; outside of dynare) to solve for your steady-state. If it did not solve within dynare, the solver should continue to iterate and not solve within matlab either. What you can then attempt is to “flip” one of the endogenous steady-state variables in your system for a parameter (i.e., turn the variable into a “parameter” and the “parameter” into a variable). If done with caution, the solver should tell you which values of the parameter are consistent with the system and thus you can avoid having to program a grid search (pain in the rear).

Here is an example:

VARIABLES: X, Y,
PARAMETERS: a, b

Equation one in your steady-state system would then be an equation with X, Y as endogenous steady-state values, and a, b will be your parameters.
Equation two in your steady-state system would then be the same.

Now, let’s say you are fairly certain about one of your parameters, b. But you are unsure of a. What you can do is “trick” the system by altering it to

Equation one in your steady-state system would then be an equation with X, b as endogenous steady-state values, and Y, a will be your parameters.
Equation two in your steady-state system would then be the same.

NOTE: you should be fairly certain about the target value for Y if you do the above. Look up the documentation on non-linear system solvers in matlab and the f-solve function in particular.