We provide non-asymptotic $L^1$ bounds to the normal for four well-known models in statistical physics and particle systems in $\mathbb{Z}^d$; the ferromagnetic nearest-neighbor Ising model, the supercritical bond percolation model, the voter model and the contact process. In the Ising model, we obtain an $L^1$ distance bound between the total magnetization and the normal distribution at any temperature when the magnetic moment parameter is nonzero, and when the inverse temperature is below critical and the magnetic moment parameter is zero. In the percolation model we obtain such a bound for the total number of points in a finite region belonging to an infinite cluster in dimensions $d \ge 2$, in the voter model for the occupation time of the origin in dimensions $d \ge 7$, and for finite time integrals of non-constant increasing cylindrical functions evaluated on the one dimensional supercritical contact process started in its unique invariant distribution. The tool developed for these purposes is a version of Stein's method adapted to positively associated random variables. In one dimension, letting $\boldsymbol{\xi}=(\xi_1,\ldots,\xi_m)$ be a positively associated mean zero random vector with components that obey the bound $|\xi_i| \le B, i=1,\ldots,m$, and whose sum $W = \sum_{i=1}^m \xi_i$ has variance 1, it holds that $$ d_1 \left(\mathcal{L}(W),\mathcal{L}(Z) \right) \leq 5B + \sqrt{\frac{8}{\pi}}\sum_{i \neq j} \mathbb{E}[\xi_i \xi_j] $$ where $Z$ has the standard normal distribution and $d_1(\cdot,\cdot)$ is the $L^1$ metric. Our methods apply in the multidimensional case with the $L^1$ metric replaced by a smooth function metric.