We study isolated and binary neutron stars in dynamical Chern-Simons gravity. This theory modifies the Einstein-Hilbert action through the introduction of a dynamical scalar field coupled to the Pontryagin density. We here treat this theory as an effective model, working to leading order in the Chern-Simons coupling. We first construct isolated neutron star solutions in the slow-rotation expansion to quadratic order in spin. We find that isolated neutron stars acquire a scalar dipole charge that corrects its spin angular momentum to linear order in spin and corrects its mass and quadrupole moment to quadratic order in spin, as measured by an observer at spatial infinity. We then consider neutron stars binaries that are widely separated and solve for their orbital evolution in this modified theory. We find that the evolution of post-Keplerian parameters is modified, with the rate of periastron advance being the dominant correction at first post-Newtonian order. We conclude by applying these results to observed pulsars with the aim to place constraints on dynamical Chern-Simons gravity. We find that the modifications to the observed mass are degenerate with the neutron star equation of state, which prevents us from testing the theory with the inferred mass of the millisecond pulsar J1614-2230. We also find that the corrections to the post-Keplerian parameters are too small to be observable today even with data from the double binary pulsar J0737-3039. Our results suggest that pulsar observations are not currently capable of constraining dynamical Chern-Simons gravity, and thus, gravitational-wave observations may be the only path to a stringent constraint of this theory in the imminent future.