We present results for matrix elements of $\Delta S=2$ four-fermion operators arising generically in models of new physics. These are needed to constrain such models using the measured values of $\varepsilon_K$ and $\Delta M_K$. We use lattice QCD with 2+1 flavors of improved staggered fermions on lattices generated by the MILC collaboration. We extrapolate to the continuum from three lattice spacings ranging down to $a\approx 0.045\;$fm. Total errors are $\sim 5-6%$, arising primarily from our use of one-loop matching between lattice and continuum operators. For two of the matrix elements, our results disagree significantly from those obtained using different fermion discretizations.