The temporal epidemic‐type aftershock sequence (ETAS) model is widely used to describe seismicity. However, only a few programs for estimation of the temporal ETAS model parameters are publicly available, and it is difficult to routinely apply some of them to observed data, due to initial value dependence. A robust temporal ETAS estimation program is required to routinely process earthquake catalogs or to perform Monte Carlo simulations. In this study, we developed a new Fortran program, etas_solve, which is based on the steepest descent method and calculates exact gradient and Hessian using the automatic differentiation technique. We performed numerical tests for a real earthquake catalog and for synthetic catalogs. Through the numerical tests, we found that initial value dependence of etas_solve is several orders of magnitude smaller than that of statistical analysis of point processes (SAPP), which is an R wrapper of statistical analysis of seismicity. We also found that standard error estimated from a Hessian at a maximum‐likelihood estimate (MLE) is useful to quantify the uncertainty of the MLE. Computation time of etas_solve was comparable to that of SAPP on a single‐core machine and faster on a multicore machine (around eightfold speed up on a 16‐core machine). In addition, etas_solve provides rich features that may be useful in real‐world applications, such as options to use fixed parameters and an auxiliary window in time and magnitude.