The discretized Laplacian is an N x N+2 stencil. To get an NxN stencil, a truncation has to occur which essentially embeds the boundary conditions into the stencil, making it an affine operator if non-zero BC, and a linear NxN if Neumann or Dirichlet with 0. You can then write out the boundary values using the condition. For example, if Dirichlet it’s trivial, since [bc(L);u;bc(R)] would be the N+2 points. Likewise you can do the same for Neumann, where a polynomial interpolation gives you a linear equation for the value at the boundary. Neumann0 with the simplest discretization of u’ implies that [u[1];u;u[end] has to hold.
So if you think of this extrapolation operator Q:N → Nx2 as the boundary conditions just written in the inverted form (pseudo-inverse in the space of functions that satisfy the BCs), then the operation is actually AQu. Since it’s much easier for people to think square, we give Δ = AQ on the interior, at least in the case of Neumann0 and Dirichlet0, and that defines the ODE on the interior while the boundary conditions are completely determined by these values and the chosen discretization (assuming linear BCs, like Robin).
[When non-zero, we don’t currently add the affine term which is the issue we need to address.]