This study presents a mathematical framework for calculating the diffusion of Brownian particles in generalized shear flows. By solving the Langevin equationsusing stochastic instead of classical calculus, we propose a mathematical formulation that resolves the particle mean-square displacement (MSD) at all timescales for any two-dimensional parallel shear flow described by a polynomial velocity profile. We show that at long timescales, the polynomial order of time of the particle MSD is n+2, where n is the polynomial order of the transverse coordinate of the velocity profile. We generalize the method to resolve particle diffusion in any polynomial shear flow at all timescales, including the order of particle relaxation timescale, which is unresolved in current theories. Particle diffusion at all timescales is then studied for the cases of Couette and plane Poiseuille flows and a polynomial approximation of a hyperbolic tangent flow while neglecting the boundary effects. We observe three main stages of particle diffusion along the timeline for Couette and plane Poiseuille flows and four main stages for hyperbolic tangent flow. The particle MSD is distinctly different across these stages due to different dominant physical mechanisms. Thus, higher temporal and spatial resolution for diffusion processes in shear flows may be realized, suggesting a more accurate analytical approach for the diffusion of Brownian particles.