With the vast amounts of data recorded by the CWB seismic networks, a fast and precise algorithm for near real-time earthquake locations as well as fundamentally urgent source and hazard studies of large damaging earthquakes in Taiwan becomes increasingly important. Here we present a novel method for automatic phase identification and event detection, modified from a Python-based PhasePApy package to be more computationally efficient and suitable for the high-rate seismicity and large-scale dense network in Taiwan. The performance efficiency is enhanced by substituting the Python code with a Fortran subroutine to calculate the characteristic function for phase picker. It is further improved by amalgamating the picks close enough in time but on different components, associating pairs of the amalgamated picks with the origin times of candidate events through an empirical linear relation between the P and S-P times reported in the CWB phase catalog, and isolating clustered candidate events and corresponding phase picks to determine the locations of real earthquakes. The algorithm is applied to detecting in almost real time over two thousands of events that occurred within four days accompanying the 4 February 2018 ML 5.89 foreshock and 6 February 2018 ML 6.25 Hualien mainshock, which far surpasses the capability of experienced human analysts. The temporal and spatial distributions of the detected aftershocks provide the timely and first-hand information to estimate the aftershock decay rate and verify the orientation of the fault plane and rupture extent of the foreshock and disastrous mainshock as a guide for fast risk assessment and source characteristics analysis.