How IRFs are computed for the model like this?

For example, I just have 2 variables, a, b, in the model and one exogenous shock e_m, which is iid standard normal.

a=0.6a(-1)+0.4a(+1)+b;

b=0.88*b(-1)+e_m;

it is easy to compute irf for variable “b”, but how irf is computed for variable “a” since there are lead and lag in the first equations?

Dynare computes the model solution to the system of equations you enter. This solution has a recursive state-space form from which IRFs are easily computed.