If the ARMA process is causal there is a general formula that provides the autocovariance coefficients.

Consider the causal $\text{ARMA}(p,q)$ process
$$
y_t = \sum_{i = 1}^p \phi_i y_{t-1} + \sum_{j = 1}^q \theta_j \epsilon_{t - j} + \epsilon_t,
$$
where $\epsilon_t$ is a white noise with mean zero and variance $\sigma_\epsilon^2$.
By the causality property, the process can be written as
$$
y_t = \sum_{j = 0}^\infty \psi_j \epsilon_{t - j},
$$
where $\psi_j$ denotes the $\psi$-weights.

The general homogeneous equation for the autocovariance coefficients of a causal $\text{ARMA}(p,q)$ process is
$$
\gamma (k) - \phi_1 \gamma (k-1) - \cdots - \phi_p \gamma (k-p) = 0,
\quad k \geq \max (p, q+1),
$$
with initial conditions
$$
\gamma (k) - \sum_{j = 1}^p \phi_j \gamma (k-j)
  = \sigma_\epsilon^2 \sum_{j = k}^q \theta_j \psi_{j - k},
\quad 0 \leq k < \max (p, q+1).
$$