T = NaN(sum(M_.dynamic_tmp_nbr(1:5)));
T = feval([M_.fname, '.dynamic_g4_tt'], oo_.dr.ys(l), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1);
g1 = feval([M_.fname, '.dynamic_g1'], T, oo_.dr.ys(l), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1, false);
g2 = feval([M_.fname, '.dynamic_g2'], T, oo_.dr.ys(l), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1, false);
g3 = feval([M_.fname, '.dynamic_g3'], T, oo_.dr.ys(l), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1, false);
g4 = feval([M_.fname, '.dynamic_g4'], T, oo_.dr.ys(l), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1, false);
The first two lines compute the temporary terms up to 4th order. The last four lines computes the matrices that you're interested in. The last argument equal to "false" means that temporary terms should not be recomputed.