Let me first say that I think Tobias Kienzler has done a great job of discussing the intuition behind your question in going from finite to infinite dimensions.
I'll, instead, attempt to address the mathematical content of Jackson's statements. My basic claim will be that
Whether you are working in finite or infinite dimension, writing the Schrodinger equation in a specific basis only involves making definitions.
To see this clearly without having to worry about possible mathematical subtleties, let's first consider
Finite dimension
In this case, we can be certain that there exists an orthnormal basis $\{|n\rangle\}_{n=1, \dots N}$ for the Hilbert space $\mathcal H$. Now for any state $|\psi(t)\rangle$ we define the so-called matrix elements of the state and Hamiltonian as follows:
\begin{align}
\psi_n(t) = \langle n|\psi(t)\rangle, \qquad H_{mn} = \langle m|H|n\rangle
\end{align}
Now take the inner product of both sides of the Schrodinger equation with $\langle m|$, and use linearity of the inner product and derivative to write
\begin{align}
\langle m|\frac{d}{dt}|\psi(t)\rangle=\frac{d}{dt}\langle m|\psi(t)\rangle=\frac{d\psi_m}{dt}(t)
\end{align}
The fact that our basis is orthonormal tells us that we have the resolution of the indentity
\begin{align}
I = \sum_{m=1}^N|m\rangle\langle m|
\end{align}
So that after taking the inner product with $\langle m|$, the write hand side of Schrodinger's equation can be written as follows:
\begin{align}
\langle m|H|\psi(t)\rangle
= \sum_{m=1}^N\langle n|H|m\rangle\langle m|\psi(t)\rangle
= \sum_{m=1}^N H_{nm}\psi_m(t)
\end{align}
Equating putting this all together gives the Schrodinger equation in the $\{|n\rangle\}$ basis;
\begin{align}
\frac{d\psi_n}{dt}(t) = \sum_{m=1}^NH_{nm}\psi_m(t)
\end{align}
Infinite dimension
With an infinite number of dimensions, we can choose to write the Schrodinger equation either in a discrete (countable) basis for the Hilbert space $\mathcal H$, which always exists by the way since quantum mechanical Hilbert spaces all possess a countable, orthonormal basis, or we can choose a continuous "basis" like the position "basis" in which to write the equation. I put basis in quotes here because the position space wavefunctions are not actually elements of the Hilbert space since they are not square-integrable functions.
In the case of a countable orthonormal basis, the computation performed above for writing the Schodinger equation in a basis follows through in precisely the same way with the replacement of $N$ with $\infty$ everywhere.
In the case of the "basis" $\{|x\rangle\rangle_{x\in\mathbb R}$, the computation above carries through almost in the exact same way (as your question essentially shows), except the definitions we made in the beginning change slightly. In particular, we define functions $\psi:\mathbb R^2\to\mathbb C$ and $h:\mathbb R^2\to\mathbb C$ by
\begin{align}
\psi(x,t) = \langle x|\psi(t)\rangle, \qquad h(x,x') = \langle x|H|x'\rangle
\end{align}
Then the position space representation of the Schrodinger equation follows by taking the inner product of both sides of the equation with $\langle x|$ and using the resolution of the identity
\begin{align}
I = \int_{-\infty}^\infty dx'\, |x'\rangle\langle x'|
\end{align}
The only real mathematical subtleties you have to worry about in this case are exactly what sorts of objects the symbols $|x\rangle$ represent (since they are not in the Hilbert space) and in what sense one can write a resolution of the identity for such objects. But once you have taken care of these issues, the conversion of the Schrodinger equation into its expression in a particular "representation" is just a matter of making the appropriate definitions.
There is nothing to "derive" for the kinetic energy operator. By definition, classical kinetic energy is $\frac{p^2}{2m}$, and so $\hat{E}_\text{kin} = \frac{\hat{p}^2}{2m}$ quantumly. It's not exactly clear why you think this doesn't make sense mathematically, but it does: In words, it says "apply the momentum operator twice, then divide the result by $2m$".
Note that $\langle p^2\rangle \neq \langle p\rangle^2$, the difference is precisely what the standard deviation is defined as and what is usually called the "uncertainty" $\Delta p$ in most physics texts.
Best Answer
The mathematics of QM seems ad-hoc when one sees it in physics courses. In von Neumann's approach it becomes clear that a good mathematical candidate for the description of quantum phenomena are operator algebras, specifically C*-algebras and von Neumann algebras.
Their representation theory then leads to the usual Dirac/Schrödinger/Heisenberg picture of quantum mechanics (it comes from the fact that, thanks to von Neumann's theorem, any system with finitely many degrees of freedom has only one irreducible representation, namely Schrödinger's, up to unitary equivalence). States are, in this picture, positive and normalised linear maps from the C*-algebra of observables to the field $\mathbb C$. Using the GNS construction, every pure state gives an irreducible representation where such state is represented by a cyclic vector.
Now the projective Hilbert space of this representation is in one-to-one correspondence with all the other states in the same superselection sector, and this justifies rigorously the use of "unit energy" vectors. Thus one can do superposition of states through normalised linear combinations, which shows the existence of the phenomenon of "quantum interference", which is absent in classical mechanics because of the commutativity of the C*-algebra of observables.