Preliminary remark. As mentioned in the comments, I find the notion "resolvent formalism", as well as the description in the Wikipedia article, rather misleading - resolvents are not somekind of formalism, and they are certainly not a mere "technique for applying complex analysis to spectral theory" (as claimed in the Wikipedia article). On the contrary, resolvents lie at the very foundation of spectral theory and are an integral part of it. Even the most basic results in spectral theory on Banach spaces - say e.g. that the spectrum of an operator is always closed - rely on the properties of the resolvent.
Answer to the question. I think the best way to understand why resolvents are so useful and important in spectral theory is to discuss several situations where resolvents occur. I do this in the following list.
Notational remark. I am going to use the sign convention $R(\lambda,A) := (\lambda-A)^{-1}$, in contrast to the Wikipedia article that defines the resolvent at $\lambda$ as $(A-\lambda)^{-1}$. The convention used in the Wikipedia article is, for instance, used in Kato's book; but aside from this book I am under the impression that the other convention is much more common in the literature.
Setting. Let $X$ be a complex Banach space and let $A: X \supseteq D(A) \to X$ be a closed linear operator. I will call $A$ bounded if $D(A) = X$ (in which case $A$ is automatically a bounded linear operator $X \to X$).
Resolvents as solution operators. We often want to solve linear equations, i.e. equations of the type
$$
(*) \qquad -Ax = y.
$$
If $0$ is not in the spectrum of $A$, then this equation has a unique solution which is given by $x = R(0,A)y$; so the resolvent of $A$ at the point $0$ is the solution operator to the equation $(*)$.
In many cases, we want to study more general equations than $(*)$, for instance the equation
$$
(**) \qquad (\lambda - A)x = y
$$
for a scalar $\lambda$. There are various reasons why we could be interested in equations of the type $(**)$; for instance:
It is a very simple generalisation of the equation $(*)$, thus raising $(*)$ to a slightly higher level of complexity.
Equations of the type $(**)$ occur in various situations; for instance, if you consider an inhomogenious linear ODE of first order, you will end up with an equation of the type $(**)$ where $\lambda = 1$.
If $0$ is in the spectrum of $A$, then the equation $(*)$ cannot be uniquely solved, in general. So we might be interested in slightly perturbing $(*)$ with a simple perturbation and solve the perturbed equation, instead. (This would lead to $(**)$ with small values of $\lambda$).
If $A$ is a differential operator, then people who work in PDE might often call the solution operator of $(**)$ (or rather, its integral kernel) the Green function of the PDE - so from a PDE perspective resolvent is an operator theoretic notion for Green function.
Locating the spectrum. The spectrum $\sigma(A)$ of $A$ consists of those complex scalars $\lambda$ for which $\lambda - A: D(A) \to X$ is not bijective.
Determining the spectrum of an operator is often an important and very often a very difficult task. But if we cannot compute the spectrum in detail we would at least like to know some information about its location.
Now an obvious way to show that a number $\lambda$ is not in the spectrum is to write down a formula for the resolvent $R(\lambda,A)$ and to show that it is indeed inverse to $\lambda - A$. This is less absurd than it may seem at first glance; here are two classical situations where one does precisely this:
Assume that $A$ is bounded and that $|\lambda| > \|A\|$. Then it follows that the so-called Neumann series
$$
\sum_{k=0}^\infty \frac{A^k}{\lambda^{k+1}}
$$
converges absolutely with respect to the operator norm, and that this operator is inverse to $\lambda - A$ (i.e. the Neumann series equals the resolvent $R(\lambda,A)$). This proves that no spectral value of $A$ has modulus strictly larger then $\|A\|$; in particlar, the spectrum of $A$ is bounded.
If $\lambda$ is not in the spectrum of $A$, then one can use a similar series expansion argument to show that numbers which are sufficiently close to $\lambda$ are not in the spectrum of $A$, either. Hence, the resolvent set of $A$ is open and thus the spectrum is closed.
Another way to use resolvents in order to obtain information about the location of the spectrum is as follows:
It is not difficult to show that the norm of the resolvent $R(\lambda,A)$ explodes as $\lambda$ approaches to spectrum of $A$. This can be useful to show that certain numbers are spectral values.
For instance, such an argument is typically used to show that the spectral radius of a positive operator on an ordered Banch space is - under certain assumptions on the space - always a spectral value. In finite dimensions, this can be used in a proof of the so-called Perron-Frobenius theorem (but there also exist other finite-dimensional proofs of this fact which do not use the resolvent).
Resolvents and functional calculi This is essentially the point of view where we consider the resolvent $R(\lambda,A)$ as the result that we obtain when we "substitute" the operator $A$ into the function $a \mapsto \frac{1}{\lambda - a}$. But why shoud we be interested in doing this?
Well, by a functional calculus we mean a method to substitute our operator $A$ into a "sufficiently nice" function $f: \mathbb{C} \supseteq \Omega \to \mathbb{C}$. For instance, if $f$ is a polynomial, this is straightforward. If $f$ is entire and $A$ is bounded, this can also be done in a straightforward way by using the Taylor series expansion of $f$.
For more general holomorphic functions, it is illuminating to recall the Cauchy integral formula: if $f: \mathbb{C} \supseteq \Omega \to \mathbb{C}$ is holomorphic, $a$ is a point in $\Omega$ and $\gamma$ is a cycle in $\Omega$ that encloses $a$ precisely once, but does not enclose any point outside of $\Omega$, then
$$
(C) \qquad f(a) = \frac{1}{2\pi i} \int_\gamma f(\lambda) \frac{1}{\lambda-a} \, d\lambda.
$$
Now if we "formally" replace $a$ with $A$ in this formula, this suggests that we define the operator $f(A)$ as
$$
(C') \qquad f(A) := \frac{1}{2\pi i} \int_\gamma f(\lambda) R(\lambda,A) \, d\lambda.
$$
In the Cauchy integral formula we required that the path $\gamma$ encloses the point $a$ - which is the only point where the function $\lambda \mapsto \frac{1}{\lambda - a}$ is not defined - exactly once. Similary, we shall assume for $(C')$ that the path $\gamma$ encloses each point $a$ where the function $\lambda \mapsto R(\lambda,A)$ is not defined - i.e., each spectral value of $A$ - exactly once.
Of course, as discussed here this is just heuristics. But it turns out that this definition of the holomorphic functional calculus works very well and has good properties.
Spectral projections. This can be seen as an application of the holomorphic functional calculus. For the sake of simplicity, assume that $A$ is bounded.
If $A$ is a finite-dimensional matrix, the underlying space $\mathbb{C}^d$ can be decomposed into generalized eigenspaces of $A$; this yields the well-known Jordan normal form of $A$. But for operators on infinite-dimensional spaces, we do not have such a normal form, in general. However, we can sometimes still decompose our operator $A$ into "simpler" operators:
To this end, assume that the spectrum of $A$ is not connected, i.e. it can be written as $\sigma(A) = \sigma_1 \cup \sigma_2$ for two non-empty and compact disjoint sets $\sigma_1$ and $\sigma_2$. Then we can find two open sets $\Omega_1$ und $\Omega_2$ which do not intersect and such that $\Omega_1$ contains $\sigma_1$ and $\Omega_2$ contains $\sigma_2$. Since our two open sets do not intersect, the indicator function $\mathbb{1}_{\Omega_1}$ is a holomorphic mapping on the union $\Omega := \Omega_1 \cup \Omega_2$.
Hence, we can use the functional calculus to define the operator $P := \mathbb{1}_{\Omega_1}(A)$, and similarly, $Q := \mathbb{1}_{\Omega_2}(A)$.
By using standard properties of the functional calculus one gets
$$
P^2 = \mathbb{1}_{\Omega_1}(A) \mathbb{1}_{\Omega_1}(A) = \Big(\mathbb{1}_{\Omega_1} \mathbb{1}_{\Omega_1}\Big)(A) = \mathbb{1}_{\Omega_1}(A) = P,
$$
and similarly, $Q^2 = Q$. Hence, $P$ and $Q$ are projections. Moreover, they sum up to
$$
P + Q = \Big(\mathbb{1}_{\Omega_1} + \mathbb{1}_{\Omega_2}\Big)(A) = \mathbb{1}_{\Omega}(A) = \operatorname{id}_X,
$$
so $P$ and $Q$ are even complementary projections. Both of them commute with $A$, so the decomposition
$$
X = PX \oplus QX
$$
reduces $A$. The projections $P$ and $Q$ are the so-called spectral projections of $A$ associated with $\sigma_1$ and $\sigma_2$. Their connection to resolvents is given by the functional calculus which we used to define $P$ and $Q$: for instance,
$$
P = \frac{1}{2\pi i} \int_\gamma \mathbb{1}_{\Omega_1}(\lambda) R(\lambda,A) \, d\lambda
$$
for an appropriately chosen path $\gamma$.
Evolutions equations: the Hille-Yosida theorem. Let $A$ be densely defined, by which we mean that the domain $D(A)$ is dense in $X$. Consider initial value problem
$$
(IVP) \quad
\begin{cases}
\dot x(t) & = Ax(t) \quad \text{for } t \ge 0, \\
x(0) & = x_0
\end{cases}
$$
for an intial vector $x_0 \in X$. The Hille-Yosida theorem say that the evolution problem (IVP) is well-posed (in a certain sense) for each $x_0 \in X$ if and only if the resolvent $R(\lambda,A)$ exists for all $\lambda$ of sufficiently large real part and its powers satisfy a certain growth estimate.
In other words: In order to show well-posedness of the evolution problem (IVP) we need to study the stationary problem $(\lambda - A)x = y$ and estimate its solutions.
So the resolvent occurs when we study well-posedness of evolution equations.
The resolvent as Laplace transform. This is related to the previous point. If $A$ is the generator of a $C_0$-semigroup $(T(t))_{t \ge 0}$ on $X$ (which means precisely that the evolution problem (IVP) is well-posed for each initial value), then for all $\lambda$ with sufficiently large real part the formula
$$
R(\lambda,A)x = \int_0^\infty e^{-\lambda t} T(t)x \, dt
$$
holds for each $x \in X$. This mean that the resolvent of $A$ is the Laplace transform of the semigroup. As a consequence of this observation, one can use so-called Tauberian theorems to obtain information about the long-term behaviour of the semigroup from information about the resolvent $R(\lambda,A)$.
An algebraich point of view: The resolvent identity. The resolvent of $A$ satisfies the so-called resolvent identity
$$
(RI) \quad (\mu - \lambda) R(\lambda,A) R(\mu,A) = R(\mu,A) - R(\lambda,A)
$$
for all scalars $\lambda, \mu$ outside the spectrum of $A$.
Now one can go one step further and call a mapping $R$ from a non-empty open subset of $\mathbb{C}$ into the space of bounded linear operators on $X$ a pseudo-resolvent if it satisfies the resolvent identity (RI).
If the operator $R(\lambda): X \to X$ is injective for one number $\lambda$ (in which case it is automatically injective for all further numbers $\lambda$, too), then one can show that there exists an operator $A$ such that $R(\lambda) = R(\lambda,A)$ for all $\lambda$.
If none of the operators $R(\lambda)$ is injective, one can still interpret $R$ as a resolvent - but of a so-called multi-valued operator rather than of a usual operator.
I do not discuss any specific applications of pseudo-resolvents here (there are applications of them!), but I thought it might be a good idea to mention them here in order to give yet one more perspective that one can take when looking at resolvents.
