I have written a function that can implement the Matrix Numerov method in a similar manner to Pillai, et al at UW-Madison, except that this script can quickly be generalized to other potential functions and run in two lines of code. I have tested it on the V potential which can also be solved by Airy functions. Thus, you can accurately approximate eigenvalues and eigenvectors for a variety of potentials.

First, a Derivation:

The Numerov algorithm is for solving differential equations of the form $$y''(x)=-g(x)*y(x) + s(x)$$ Invoke a Taylor expansion of y(x) about $x_0:$ $$y(x)= y(x_0) +(x-x_0)\; y'(x_0) +\frac{(x-x_0)^2 }{2!} y''(x_0) +\frac{(x-x_0)^3}{3!} y'''(x_0) +\frac{(x-x_0)^4}{4!} y^{iv}(x_0)+...$$ Rewrite, letting $h$ equal the distance from point $x_0$ to $x$ such that $y(x)=y(x0+h)$: $$y(x_0+h) =y(x_0) +h\;y'(x_0) +\frac{h^2}{2!} y''(x_0) +\frac{h^3}{3!}\;y'''(x_0) +\frac{h^4}{4!} y^{iv}(x_0)+...$$ Define a grid of x values and let $(x0+h)=x_{n+1}$ and $x0=x_n$. Then a step forward (toward positive x) is: $$y_{n+1} =y_{n}+ h\;y'(x_{n} ) +\frac{h^2}{2!}\;y''(x_{n}) +\frac{h^3}{3!}\;y'''(x_{n})+ \frac{h^4}{4!} y^{iv}(x_n)+...$$ That was a step forwards. For a step backwards we replace $h$ with $-h$: $$y_{n-1} =y_{n}- h\;y'(x_{n} ) +\frac{h^2}{2!}\;y''(x_{n}) -\frac{h^3}{3!}\;y'''(x_{n}) +\frac{h^4}{4!} y^{iv}(x_n)+...$$ Add the two equations together and you get: $$y_{n+1} +y_{n-1} =2 y_{n} +h^2 y''(x_{n}) +\frac{h^{4}}{12} y^{iv}(x_n) +O(h^6)$$ Rearrange: $$y_{n+1} -2y_n +y_{n-1} =\mathbf {h^2 y''(x_n)} +\mathbf{\frac{h^4}{12} y^{iv}(x_n)}+O(h^6)$$ We will largely ignore $O(h^6)$ and substitute for the bold $h^2 y''(x_n)$ and $\frac{h^4}{12} y^{iv}(x_n)$ above. For $\mathbf{h^2 y''(x_n)}$ we use the original differential equation: $y''(x_n) =-g_n y_n +s_n$.

For the bold $\mathbf{\frac{h^4}{12} y^{iv}(x_n)}$, we figure that:

$$\frac{h^4}{12}\frac{d^2}{dx^2} (y''(x_n)) =\frac{h^4}{12}\frac{d^2}{dx^2} (-g_n y_n +s_n)$$ Now we just learned (rearranging the prior $y_{n+1} -2y_n +y_{n-1}$ equation again), that for generic $y(x)$, $$y''(x_n) =\frac{y_{n+1} -2y_n +y_{n-1}-O(h^4)}{h^2}$$ Use that idea twice to express the pieces of $\frac{d^2}{dx^2} (-g_n y_n +s_n)$ and you get: $$\frac{h^4}{12}\frac{d^2}{dx^2} (-g_n y_n +s_n) =\frac{h^4}{12} \left[ \frac{-(g_{n+1} y_{n+1} -2g_n y_n +g_{n-1} y_{n-1})}{h^2} +\frac{s_{n+1} -2s_n +s_{n-1}} {h^2} \right]$$ ...which is then equal to: $$\frac{h^2}{12} [-(g_{n+1} y_{n+1} -2g_{n} y_{n} +g_{n-1} y_{n-1}) +s_{n+1} -2s_{n} +s_{n-1}]$$ Now sub that back into the equation with the boldface: $$y_{n+1} -2y_n +y_{n-1} =h^2(-g_n y_n +s_n) +\frac{h^2}{12} [-g_{n+1} y_{n+1} +2g_{n} y_{n} -g_{n-1} y_{n-1} +s_{n+1} -2s_{n} +s_{n-1}]$$ Divide both sides by $h^2$: $$\frac{y_{n+1} -2y_n +y_{n-1}}{h^2} =(-g_n y_n +s_n) +\frac{1}{12} [-g_{n+1} y_{n+1} +2g_{n} y_{n} -g_{n-1} y_{n-1} +s_{n+1} -2s_{n} +s_{n-1}]$$ ...which is equal to (combining terms): $$\frac{-1}{12}(g_{n+1} y_{n+1} +10g_n y_n +g_{n-1} y_{n-1}) +\frac{1}{12} (s_{n+1} +10s_n +s_{n-1})$$ Now, applying the regular Numerov algorithm might consist of picking two points $y_{n-1}$ and $y_n$, using the above formula to generate $y_{n+1}$ and iterating over all points $n$ in the $x$ grid. But here I will substitute into the one-dimensional time-independent Schrödinger equation and introduce the matrix method. Rearrange the Schrödinger Equation: The one-dimensional Schrödinger equation gives $$\frac{-\hbar^2}{2m}\frac{d^2\psi}{dx^2} +V(x)\psi(x)=E\psi(x)$$ Rearranging, $$\frac{d^2\psi}{dx^2}= \frac{2m}{\hbar^2}V(x)\psi(x) -\frac{2mE}{\hbar^2}\psi(x)$$ Looking familiar yet? Let $y''(x)=\frac{d^2\psi}{dx^2}$, and $-g(x)=\frac{2m}{\hbar^2}V(x)$, $y(x)=\psi(x)$, and $s(x)=\frac{-2mE}{\hbar^2}\psi(x)$. Now again define a grid of x values (how to get the grid spacing right to come later...) and let $\psi_n$ represent the value of $\psi(x)$ at the $x_n$'th point. Now we can use the result above, noting that $\hbar$ is the reduced Planck constant, but $h$ is the grid spacing $\Delta x$: $$\begin{align} \frac{y_{n+1} -2y_n +y_{n+1}}{h^2} = & \frac{-1}{12} \left[ (\frac{-2m}{\hbar^2}V_{n+1}) (\psi_{n+1} ) +10(\frac{-2m}{\hbar^2} V_n) \psi_n +(\frac{-2m}{\hbar^2}V_{n-1}) (\psi_{n-1} )\right] \\\ & +\frac{1}{12} \left[ \frac{-2mE}{\hbar^2}\psi_{n+1} +10(\frac{-2mE}{\hbar^2})\psi_n +\frac{-2mE}{\hbar^2}\psi_{n-1} \right] \end{align}$$ Multiply everything by $\frac{-\hbar^2}{2m}$ and rearrange a little: $$\frac{-\hbar^2}{2m}\left[ \frac{\psi_{n+1} -2\psi_n +\psi_{n-1}}{h^2} \right] +\frac{1}{12} \left[V_{n+1} \psi_{n+1} +10V_n \psi_n +V_{n-1} \psi_{n-1}\right] =\frac{1}{12}E \left[ \psi_{n+1} +10\psi_n +\psi_{n-1} \right]$$ And now for...

The Matrix Form

Let $\psi(x)$ be represented by a column vector... I.e. $$ \begin{pmatrix} \psi_{1} \\\ \psi_{2} \\\ \psi_{3} \\\ \vdots \end{pmatrix} $$

Here each $\psi_i$ is the value of $\psi(x)$ at the grid point $x_i$ (Many apologies if the above shows up as a row vector; I tried). Now we can express the equation above in matrix format as... $$\frac{-\hbar^2}{2m}\mathbf{A} \;\mathbf{\psi} +\mathbf{B}\;\mathbf{V}\;\mathbf{\psi} =E\;\mathbf{B}\;\mathbf{\psi}$$ Now, $\mathbf{A}$ is equal to... $\frac{1}{h^2}*$ $$ \begin{pmatrix} -2 & 1 & 0 & 0 & 0 &\cdots\\ 1 & -2 & 1 & 0 & 0 &\cdots\\ 0 & 1 & -2 & 1 & 0 &\cdots\\ 0 & 0 & 1 & -2 & 1 &\cdots\\ 0 & 0 & 0 & 1 & -2 &\cdots\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix} $$

And B is $\frac{1}{12}*$ $$ \begin{pmatrix} 10 & 1 & 0 & 0 & 0 &\cdots\\ 1 & 10 & 1 & 0 & 0 &\cdots\\ 0 & 1 & 10 & 1 & 0 &\cdots\\ 0 & 0 & 1 & 10 & 1 &\cdots\\ 0 & 0 & 0 & 1 & 10 &\cdots\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix} $$

Finally, V is equal to the diagonal matrix: $$ \begin{pmatrix} V_{1} & 0 & 0 & 0 & 0 &\cdots\\ 0 & V_{2} & 0 & 0 & 0 &\cdots\\ 0 & 0 & V_{3} & 0 & 0 &\cdots\\ 0 & 0 & 0 & V_{4} & 0 &\cdots\\ 0 & 0 & 0 & 0 & V_{5} &\cdots\\ \vdots & \vdots & \vdots & \vdots & \vdots &\ddots \end{pmatrix} $$ (Does this remind anyone of Discrete Variable Representation? Mathematica code for that to come soon...)

Left-multiplying by $\mathbf{B^{-1}}$ gives... $$\frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A\;\psi} +\mathbf{V\;\psi} =E\;\psi$$ Rearranging, $$\left( \frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A} +\mathbf{V} \right) \psi=E\;\psi$$ And voila! we have an equation in the $\mathbf{H} \psi = E \psi$ form. So $\frac{-\hbar^2}{2m}\mathbf{B^{-1}\;A}$ is your kinetic energy operator. If you're not sure the matrix form is equivalent, pick three elements of the $\psi$ vector to be $\psi_{n-1}$, $\psi_n$, and $\psi_{n+1}$, and propagate through the matrix algebra element-wise.

I will now define these matrices and operations in Mathematica...

Mathematica Code

Grid spacing should be fine enough to have five points in the narrowest half-wavelength. For the problem I was working on, I cared about wavefunctions near an energy of 8.84755*10^-16 ergs, so $p0 = \sqrt{2 m E} = 1.4327*10^-19$ --> the narrowest lobe has width $\frac{1}{2}*\frac{h}{p0}=2.31244*10^-8$ (this time h = Planck's constant in ergsec). The grid spacing should be $\frac{1}{5}$ of that width, or 4.6248810^-9... I believe this problem was using inches for some unfathomable reason.

To use this function, first define your potential function as <some name> = <some function>, i.e.:

v1input[x_] = a Abs[x]

(Here a was a constant set equal to 5.368244090805358*10^-9, the slope of the V-shaped potential). Then, decide your minimum and maximum values of x (advice for the V-shaped potential: go a couple wavelengths beyond the wall at the energy of interest...).

Then employ the following function, where expr_ is the name of your potential function...

MatrixNumerov[expr_, xmin_, xmax_, dx_, hbar_, m_] := 
 Module[{xvals, Vfunction, Vmatrix, A, B, KE, n1, H, eivals, eivecs},
  xvals = Table[i, {i, xmin, xmax, dx}];
  Vfunction = expr[xvals];
  Vmatrix = DiagonalMatrix[Vfunction];
  n1 = Length[xvals];
  B = (DiagonalMatrix[1 + 0 Range[n1 - Abs[-1]], -1] + 
      10 DiagonalMatrix[1 + 0 Range[n1 - Abs[0]], 0] + 
      DiagonalMatrix[1 + 0 Range[n1 - Abs[1]], 1])/12.; 
  A = (DiagonalMatrix[1 + 0 Range[n1 - Abs[-1]], -1] - 
      2 DiagonalMatrix[1 + 0 Range[n1 - Abs[0]], 0] + 
      DiagonalMatrix[1 + 0 Range[n1 - Abs[1]], 1])/((dx)^2); 
  KE = (-hbar^2/(2 m))*Inverse[B].A;
  H = KE + Vmatrix;
  {eivals, eivecs} = Eigensystem[H];
  {eivals, eivecs}]

Here's an example... With x running from -2.573102897338254`E-7 to 2.573102897338254E-7 in increments of 4.624881629283604'E-9,

the code:

MatrixNumerov[v1input, -2.573102897338254`*^-7, 
2.573102897338254`*^-7, 4.624881629283604`*^-9, 1.0545718`*^-27, 
1.1600000000000001`*^-23]

...returns a list of eigenvalues and a list of eigenvectors, just like a regular Eigensystem[...] command, with the smallest ten eigenvalues equal to $8.84669*10^{-16}, 8.21013*10^{-16}, 7.55796*10^{-16}, 6.86411*10^{-16}, 6.14806*10^{-16}, 5.36843*10^{-16}, 4.5527*10^{-16}, 3.61808*10^{-16}, 2.60396*10^{-16}, 1.13641*10^{-16}$ . Those ten eigen-energies solved by the Airy method were: $8.84755*10^{-16}, 8.21055*10^{-16}, 7.5585*10^{-16}, 6.8642*10^{-16}, 6.14836*10^{-16}, 5.36824*10^{-16}, 4.55283*10^{-16}, 3.61758*10^{-16}, 2.604*10^{-16}, 1.13465*10^{-16}$ So you can see, the approximation is darn close. This has been adapted from my quantum mechanics midterm (MIT course 5.73 taught by Bob Field Fall 2018) - apologies for any errors/omissions, and I'm afraid I'll have to leave a discussion of how to get the energies by the Airy functions for V-shaped potential with slope = a (above) for another time.