The number of data points for spline interpolation is $n+1$. If you want to interpolate the function $f$ on the grid $x_0,x_1,\ldots,x_n$, the only data you use is the value of $f$ on the grid points, that is, you require the spline function $s$ to satisfy
$$
s(x_i) = f(x_i), \qquad \text{for } i=0,1,\ldots,n.
$$
Additionally, the spline function is required to be $C^2$ (twice continuously differentiable). That is, if $s_i$ denotes the restriction of $s$ to the interval $[x_i,x_{i+1}]$, then
$$
s_{i-1}(x_i) = s_i(x_i) \quad\text{and}\quad
s'_{i-1}(x_i) = s'_i(x_i) \quad\text{and}\quad
s''_{i-1}(x_i) = s''_i(x_i) \qquad\text{for } i=1,2,\ldots,n-1.
$$
This gives $3(n-1)$ additional conditions, which together with the $n+1$ interpolation conditions yield the total of $4n-2$ quoted by the original poster.
However, the only data used in the $4n-2$ conditions is $f(x_i)$ for $i=0,1,\ldots,n$ and this is $n+1$ pieces of data. The continuity conditions do not use any data from the function $f$.
Definitions
L
$L$ is not the length of the interval, but a Lipschitz constant: some
number such that $$|f(x,y_1) - f(x,y_2)| \le L|y_1 - y_2|$$
satisfied for all $y_1,y_2,x$ (in some interval of interest).
M
Similarly, $M$ is any number such that the assumption $|y''(x)| \le M$
is satisfied. The error analysis is also explained in this article
here (reference taken from the Wikipedia article on the Euler
method).
Example
Suppose the equation $y' = -\frac12 y$ with initial condition $y(0) = 1$ with the Euler method on the interval $x \in [0,1]$. We actually know the exact solution in this case, $y(x) = \exp(\frac12x)$.
Finding L
We note that $f(x,y) = \frac12 y$ in the example, so the equation
$$|f(x,y_1) - f(x,y_2)| \le L|y_1 - y_2|$$
becomes
$$|\frac12 y_1 - \frac12 y_2| \le L |y_1 - y_2|.$$
This equation is satisfied for $L = \frac12$, so this is the value of $L$ we will take (we can also have taken a bigger value for $L$, like $L=1$, in which case we'll get a result which is also valid but not as sharp).
Finding M
We need it to satisfy $|y''(x)| \le M$. Since we know the exact solution, $y(x) = \exp(\frac12x)$, we know that $y''(x) = \frac14\exp(\frac12x)$ and so $|y''(x)| \le \frac14$ on the interval $x \in [0,1]$. However, for most equations we do not know the exact solution, so let us pretend we do not. Then what we can do is the following: differentiate the differential equation $y'(x) = \frac12 y$ to get $y''(x) = \frac12 y'$, and then substitute the differential equation in it to get
$$y''(x) = \frac12 \cdot \frac12 y = \frac14y.$$
We do know that the numerical solution given by the Euler method is a decreasing sequence for this example, so $y\le1$ and thus $|y''(x)| \le \frac14$ (as we found before).
Conclusion
So, we can take $L = \frac12$ and $M = \frac14$ in the error bound.
Best Answer
The direct error in your formula is that
Indeed, as reconstructed below, you get a correct bound with the same terms, but in a sum, not a difference.
To simplify the calculus, consider an autonomous system $y'=f(y)$, as you can always have $f_1(y)=1$ to get a time variable. The discretization error of Heun's explicit trapezoidal method is, in the form you use, given by \begin{align} E(h)&=y(x+h)-y(x)-\frac{h}2\Bigl(f(y(x))+f\bigl(y(x)+hf(y(x))\bigr)\Bigr) \\ &=\left[y(x+h)-y(x)-\frac12\bigl(f(y(x))+f(y(x+h))\bigr)\right]+\frac{h}2\left[f\bigl(y(x)+hf(y(x))\bigr)-f(y(x+h))\right]. \end{align}
For the first term use the error formula of the trapezoidal quadrature method for the function $g(x)=y_i'(x)=f_i(y(x))$. It gives $$ \frac{h}2(g(x)+g(x+h))-\int_x^{x+h}g(s)ds=-\frac{h^3}{12}g''(x+\theta_i h) $$ so that the first term has the error bound \begin{align} \left\|y(x+h)-y(x)-\frac h2(f(y(x+h))+f(y(x)))\right\| &\le\frac{h^3}{12}\sup_{\theta\in[0,1]}\|y'''(x+\theta h)\| \end{align} For the second term we get using a Lipschitz constant $L_f=\|f'\|_\infty$ \begin{align} \|f(y(x)+hf(y(x)))−f(y(x+h))\|&\le L_f\|y(x)+hy'(x)−y(x+h)\| \\&\le L_f\int_x^{x+h}(x+h-s)\|y''(s)\| \\&\le L_f \frac{h^2}2\sup_{\theta\in[0,1]}\|y''(x+\theta h)\| \end{align}
In combination this gives for the error of the explicit method \begin{align} \|E(h)\| &\le \frac{h^3}{12}·\|y'''\|_\infty+\frac{h^3}{4}·\|f'\|_\infty·\|y''\|_\infty \end{align}
This is the most optimistic outcome you could expect from the error formula of the linked article.