Let $A\in\mathbb{R}^{n\times n}$ a diagonalizable stable matrix (all eigenvalues have a negative real part) and consider the difference of its resolvent matrix evaluated at symmetric points on the complex line \begin{align} &(j\omega I -A)^{-1} - (-j\omega I -A)^{-1} \cr =& (j\omega I -A)^{-1} + (j\omega I +A)^{-1} \tag{a} \end{align} I have arrived to the following expression $$ (j\omega I -A)^{-1} + (j\omega I +A)^{-1} = -2j\omega (I\omega^2 +A^2)^{-1} \tag{b} $$ Is the following proof ok?
Proof: Because $A$ have no eigenvalues on the pure imaginary line, the expression in (a) is analytic. Now we use the following expression to evaluate matrix analytic functions \begin{equation} f(A) = \frac{1}{2j\pi}\int_{\Gamma} (zI-A)^{-1} f(z) dz \tag{c} \end{equation} where $\Gamma$ is a closed contour which encloses the spectrum of $A$. We have \begin{align} (j\omega I -A)^{-1} + (j\omega I +A)^{-1} &= \frac{1}{2j\pi}\int_{\Gamma} (zI-A)^{-1} \left(\frac{1}{j\omega - z} +\frac{1}{j\omega + z}\right) dz\cr &= \frac{1}{2j\pi}\int_{\Gamma} (zI-A)^{-1} \left(\frac{-(j\omega+z)-(j\omega-z)}{\omega^2 + z^2}\right) dz\cr &=\frac{1}{2j\pi}\int_{\Gamma} (zI-A)^{-1} \left(\frac{-2j\omega}{\omega^2 + z^2}\right) dz \cr &=-2j\omega(I\omega^2 + A^2)^{-1} \end{align} It is quite strange to manipulate matrices as scalars, but the proof seems good for me. Can somebody find another proof? Is this expression right?
Here is a plot of the error for this identity made in Matlab with a random stable matrix.

It sure can get more elementary: $$ \begin{split} (j\omega I-A)^{-1}+(j\omega I+A)^{-1} &= (j\omega I-A)^{-1}[(j\omega I+A)+(j\omega I-A)](j\omega I+A)^{-1} \\&= 2j\omega(j\omega I-A)^{-1}(j\omega I+A)^{-1} \\&= 2j\omega[(j\omega I+A)(j\omega I-A)]^{-1} \\&= 2j\omega[-\omega^2I-A^2]^{-1} \\&= -2j\omega[\omega^2I+A^2]^{-1}. \end{split} $$