The formation of small and large-scale hydrodynamic features in planetary nebulae are of great interest to scientist for their role in distributing matter to the interstellar medium. These structures contain both gas and dust but their origins and evolution remain mysterious. There has been much debate on the roles of dust and radiation in creating these hydrodynamic features; they are subject to complex two-way coupling: gas-particles, particles-radiation, and gas-radiation. The origin of these hydrodynamic perturbations, whether arising in the dust-driven winds of the AGB phase or sculpted by fluid and radiation forces in the early PNe phase, remains unknown. In this work, we broadly investigate the role of dust and radiation in the formation of small-scale hydrodynamic features known as cometary knots. Previous research studies have considered similar physics in dust-driven winds at shorter length and time scales, using either 1D simulations, or 2D simulations with a single mixed particle-gas fluid. Here we present a Eulerian-Lagrangian method for studying this problem at larger length and time scales. Simulations are performed using the FLASH code, developed at the FLASH Center at the University of Chicago. The Particle-in-Cell method was used with the two-dimensional Euler equations and solved using the directionally split piecewise-parabolic method. This method was then modified for the astrophysics regime by implementing radiation and non-continuum drag models for the particle and gas phases. The effects of a perturbed radiation field and perturbed particle spatial distribution were investigated to determine if these could be responsible for the formation of cometary knots observed in planetary nebulae.