Since this question has multiple sub-questions in edits, comments on answers, etc., and these have not been addressed, here goes.
Matched filters
Consider a finite-energy signal $s(t)$ that is the input to a (linear
time-invariant BIBO-stable) filter with impulse response $h(t)$, transfer function $H(f)$,
and produces the output
signal
$$y(\tau) = \int_{-\infty}^\infty s(\tau-t)h(t)\,\mathrm dt.\tag{1}$$
What choice of $h(t)$ will produce a maximum response at a given time
$t_0$? That is, we are looking for a filter such that the global maximum
of $y(\tau)$ occurs at $t_0$. This really is a very loosely phrased
(and really unanswerable) question because clearly the filter
with impulse response $2h(t)$ will have larger response than
the filter with impulse response $h(t)$, and so there is
no such thing as the filter that maximizes the response.
So, rather than compare apples and oranges, let us include the
constraint that we seek the filter that maximizes $y(t_0)$ subject
to the impulse response having a fixed energy, for example, subject to
$$\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt = \mathbb E
= \int_{-\infty}^\infty |s(t)|^2 \,\mathrm dt.\tag{2}$$
Here onwards, "filter" shall mean a linear time-invariant filter
whose impulse response satisfies (2).
The Cauchy-Schwarz inequality provides an answer to this question. We have
$$y(t_0) = \int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt
\leq \sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}
= \mathbb E$$
with equality occurring if $h(t) = \lambda s(t_0-t)$ with $\lambda > 0$
where from (2) we get that $\lambda = 1$, that
is, the filter with impulse response $h(t) = s(t_0-t)$ produces
the maximal response $y(t_0) = \mathbb E$ at the specified time $t_0$.
In the (non-stochastic) sense described above, this filter is
said to be
the filter matched to $s(t)$ at time $t_0$ or
the matched filter for $s(t)$ at time $t_0.$
There are several points worth noting about this result.
The output of the matched filter has a
unique global maximum value of $\mathbb E$ at $t_0$; for any other
$t$, we have $y(t) < y(t_0) = \mathbb E$.
The impulse response $s(t_0-t) = s(-(t-t_0))$
of the matched filter for time $t_0$ is just $s(t)$ "reversed in time"
and moved to the right by $t_0$.
a. If $s(t)$ has finite support, say, $[0,T]$, then the matched filter is
noncausal if $t_0 < T$.
b. The filter matched to $s(t)$ at time $t_1 > t_0$ is just the filter
matched at time $t_0$ with an additional delay of $t_1-t_0$. For this
reason, some people call the filter with impulse response $s(-t)$,
(that is, the filter matched to $s(t)$ at $t=0$) the matched filter for $s(t)$ with the
understanding that the exact time of match can be incorporated into
the discussion as and when needed. If $s(t) = 0$ for $t < 0$, then
the matched filter is noncausal. With this, we can rephrase 1. as
The matched filter for $s(t)$ produces a unique global maximum
value $y(0) = \mathbb E$ at time $t=0$. Furthermore,
$$y(t) = \int_{-\infty}^\infty s(t-\tau)s(-\tau)\,\mathrm d\tau
= \int_{-\infty}^\infty s(\tau-t)s(\tau)\,\mathrm d\tau = R_s(t)$$
is the autocorrelation function of the signal $s(t)$. It is
well-known, of course, that $R_s(t)$ is an even function of $t$
with a unique peak at the origin. Note that the output of the
filter matched at time $t_0$ is just $R_s(t-t_0)$, the autocorrelation
function delayed to peak at time $t_0$.
No filter other than the
matched filter for time $t_0$ can produce an output as large
as $\mathbb E$ at $t_0$. However, for any $t_0$,
it is possible to find filters that
have outputs that exceed $R_s(t_0)$ at $t_0$. Note that $R_s(t_0) < \mathbb E$.
The transfer function of the matched filter is $H(f)=S^*(f)$, the
complex conjugate of the spectrum of $S(f)$.
Thus, $Y(f) = \mathfrak F[y(t)]= |S(f)|^2$.
Think of this result as follows. Since $x^2 > x$ for $x > 1$ and $x^2< x$ for
$0 < x < 1$, the matched filter has low gain at those frequencies where
$S(f)$ is small, and high gain at those frequencies where $S(f)$ is large.
Thus, the matched filter is reducing the weak spectral components
and enhancing the strong spectral components in $S(f)$. (It is also
doing phase compensation to adjust all the "sinusoids" so that
they all peak at $t=0$).
-------
But what about noise and SNR and stuff like that which is what the OP
was asking about?
If the signal $s(t)$ plus additive white Gaussian noise with
two-sided power spectral density $\frac{N_0}{2}$ is processed
through a filter with impulse response $h(t)$, then the output
noise process is a zero-mean stationary Gaussian process with
autocorrelation function $\frac{N_0}{2}R_s(t)$. Thus, the
variance is
$$\sigma^2 = \frac{N_0}{2} R_s(0) = \frac{N_0}{2}\int_{-\infty}^{\infty} |h(t)|^2\,\mathrm dt.$$
It is important to note that the variance is the same regardless
of when we sample the filter output. So, what choice of $h(t)$
will maximize the SNR $y(t_0)/\sigma$ at time $t_0$? Well, from the
Cauchy-Schwarz inequality, we have
$$\text{SNR} = \frac{y(t_0)}{\sigma}
= \frac{\int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}
\leq \frac{\sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}} = \sqrt{\frac{2\mathbb E}{N_0}}$$
with equality exactly when $h(t) = s(t_0-t)$, the filter that is matched
to $s(t)$ at time $t_0$!! Note that $\sigma^2 = \mathbb EN_0/2$.
If we use this matched filter for our desired sample time, then at other
times $t_1$, the SNR will be
$y(t_1)/\sigma < y(t_0)/\sigma = \sqrt{\frac{2\mathbb E}{N_0}}$. Could
another filter give a larger SNR at time $t_1$? Sure, because $\sigma$
is the same for all filters under consideration, and we have noted above that
it is possible to have a signal output larger than $y(t_1)$ at time
$t_1$ by use of a different non-matched filter.
In short,
"does the matched filter maximize the SNR only at the sampling
instant, or everywhere?" has the answer that the SNR is maximized only
at the sampling instant $t_0$. At other times, other filters could give
a larger SNR than what the matched filter is providing at time $t_1$,
but this still smaller than the SNR $\sqrt{\frac{2\mathbb E}{N_0}}$
that the matched filter is giving you at $t_0$, and if desired,
the matched filter could be redesigned to produce its peak at time
$t_1$ instead of $t_0$.
"why not make a filter that makes a really tall skinny spike at the point of decision. Wouldn't that make the SNR even better?"
The matched filter does produce a spike of sorts at the sampling time
but it is constrained by the shape of the autocorrelation function. Any
other filter that you can devise to produce a tall skinny (time-domain)
spike is not a matched filter and so will not give you the largest possible
SNR. Note that increasing the amplitude of the filter impulse response
(or using a time-varying filter that boosts the gain at the time
of sampling) does not change the SNR since both the signal and the noise standard deviation increase proportionately.
"The I&D will basically ramp up until it gets to the sampling time and the idea is that one samples at the peak of the I&D because at that point, the SNR is a maximum."
For NRZ data and rectangular pulses, the matched filter impulse response is
also a rectangular pulse. The integrate-and-dump circuit is a correlator
whose output equals the matched filter output only at the sampling instants, and not in-between. See the figure below.

If you sample the correlator output at other times,
you get noise with smaller variance but you can't simply add up the samples
of I&D output taken at different times because the noise variables are highly correlated, and
the net variance works out to be much larger. Nor should you expect to be able
to take multiple samples from the matched filter output and combine them
in any way to get a better SNR. It doesn't work. What you have in effect
is a different filter, and you cannot do better than the (linear)
matched filter in Gaussian noise; no nonlinear processing will give
a smaller error probability than the matched fiter.