We present a quasi-three-dimensional efficient model for simulating and designing the terahertz quantum cascade laser with nonlinear axial waveguide structure, based on the finite difference beam propagation method. The traditional beam propagation method is widely used to simulate the beam profile of the passive waveguide. In order to study the active device, however, the current induced variation in the active region should also be considered in the numerical simulation model. In the model presented in this paper, the phase and the amplitude of the propagating confined field in the active waveguide are determined by a few linear and non-linear effects. The parameters relating to the linear effects, such as the intrinsic refractive index profile and the intrinsic losses of the waveguide under zero current injection, are calculated by using COMSOL-Multiphysics. While the non-linear effects, such as the modal gain and the refractive index variation induced by current injection, are considered in a rigorous way by including the rate-equation set for calculating the carrier dynamics in the active region. The parameters used in the rate-equation set are obtained by referring to the literature and fitting the experimental results of the considered terahertz lasers. By adding the current induced gain and refractive index variation, the presented beam propagation model is able to simulate many current-dependant properties of a laser, such as the output power, the gain guiding effect, and the self-focusing effect. We show in this paper that the latter two effects have influence on inner-waveguide beam profile, and the competitive balance between them determines the output beam quality. By utilizing this numerical model, the terahertz quantum cascade laser with tapered waveguide structure is simulated, and the influences of the taper angle on output power and beam quality are investigated. According to the simulation results, we find that there is an obvious increase in the output power when the taper angle is increased from 0 to 3 degree, while the increment in the output power decreases rapidly when the taper angle is further increased. Besides, we observe that for the far field the full width at half maximum of the output beam decreases sharply with increasing the taper angle. However, when the taper angle equals 8 degree, multiple lateral modes are observed, which indicates poor output beam quality of this device and poor beam coupling efficiency between this device and the power meter.Therefore, although the simulation results show that the output power of this device is higher than that of the device with 5 degree taper angle, the experiment results show that the measured output power is lower. So the taper angle is not the larger the better, but there exists an optimum value, at which the terahertz quantum cascade laser can achieve the highest effective output power.