Centered finite differences for the zero-Dirichlet boundary-value problem

The input is a real-valued continuous function $g$ on $[0,1]$ and a positive integer $N$. The procedure uses exact real arithmetic. Its admissible inputs have $N\ge 1$; an input that is not a positive integer returns invalid-mesh.

Set $h=1\mathrm{/}(N+1)$, $x_j=jh$ for $0\le j\le N+1$, and $g_j=g(x_j)$ for $1\le j\le N$. The unknowns are the $N$ interior values $U_1,\dots ,U_N$; set the boundary values $U_0=U_{N+1}=0$. The finite-difference equations are

or $A_{h}U=g_h$, where $A_h=h^{-2}\mathrm{tridiag}(-1,2,-1)$, with the known zero boundary values moved to the right-hand side (they contribute zero here).

To solve, form $b_j=h^2{g}_j$. Set $p_1=2$ and $s_1=b_1$. For $j=2,\dots ,N$, in increasing order, if $p_{j-1}=0$ return linear-solve-failure; otherwise set

If a formed pivot $p_j$ is zero, return linear-solve-failure. Otherwise set $U_N=s_N\mathrm{/}{p}_N$ and, in decreasing order for $j=N-1,\dots ,1$, set

Check the equations exactly. If any equation fails, return linear-solve-failure; otherwise return status success with the vector $U$ and the boundary values $U_0=U_{N+1}=0$.

This is an exact-arithmetic definition and makes no finite-precision accuracy claim. A numerical implementation that returns an approximate vector must report its residual separately; the assigned error theorem includes that residual as an independent term.