There is indeed such a procedure. The starting point is the pair of equations of state which express $T$ and $S$ as functions of $p$ and $V$ and so establish the model you are using---for the ideal gas, we have $T=p V, \quad S=\frac1{\gamma-1} \ln (p V^\gamma)$ (we have omitted physical constants to simplify the expressions---mathematicians tend to do that). The general case is $T=f(p,V), \quad S=g(p,V)$ for suitable functions $f$ and $g$ of two variables. $f$ and $g$ are related by the Maxwell relations which mean that $f_1g_2-f_2g_1=1$ (subscripts denote partial derivatives with respect to $p$ and $V$). Now one can express all of the basic thermodynamical quantities in terms of $p$, $V$, $f$, $g$ and the partials of the latter two. One then verifies an identity by substituting these in both sides and checking by simple algebra if it is valid. This is a simple application of the inverse function theorem and the chain rule. For those who are not comfortable with these techniques, one can write a simple Mathematica programme which does this for you and this has been carried out in the arXiv paper "A systematic approach to thermodynamical identities" (1102.1540) which, in turn, is based on the theory developed in the arXiv paper 1102.1540. The former paper even develops a method for {\it generating} thermodynamical identities, using the technique of Groebner bases.
As a small sample, one finds that the expression $C_p-C_V$ that you mention is $\dfrac f{f_1f_2}$ and from this it follows easily that this is constant for the ideal gas (for which $f(p,V)=p V$).
One sees the advantages of this method when one considers more elaborate models, e.g., the van der Waals gas ($f(p,V)=\left(p + \frac a {V^2} \right )(V-b)$, $S=\frac 1 {\gamma-1}( \ln (p+ \frac a {V^2})+\gamma \ln (V-b))$),
or the Feynman model (like the ideal gas but allowing for the fact that the adiabatic index $\gamma$ is not constant for a real gas but depends on temperature) which is analysed in the above-mentioned articles.
Added at Nick's request: $C_p$ and $C_V$ are defined as $T \frac {\partial S}{\partial T} |_V$ and $T \frac {\partial S}{\partial T} |_p$. The partial derivatives are, in the numerical code ($p\to 1$, $V \to 2$, $T \to 3$ and $S \to 4$) introduced in these articles, $(4,3,2)$ and $(4,3,1)$. This leads to the relationships $C_V= f \frac{g_1}{f_1}$ and $C_p=f\frac {g_2}{f_2}$. Taking the difference leads to the required result (since $f_1g_2-f_2g_1=1$---one of the manifestations of the Maxwell relations which can be interpreted as stating that the map from the $p,V$ plane to the $T,S$ one described by the equations of state is area preserving).
Interestingly, the difference of these two terms depends only on $f$ (i.e., temperature as a function of pressure and volume) and not on entropy, in contrast to the thermal capacities themselves. Using these formulae, we can easily compute these quantities for models which are more elborate than the ideal gas, something we have not encountered in the literature.
