Finite Difference Boundary Conditions A simple example will be the finite difference equivalent of $\partial_x^2  u(x) = b(x) $. Define $K$ and $\mathbf{u} $ as
$$K=\begin{bmatrix}
2 &-1 & 0 &0 \\
-1&2 & -1 &0 \\
0 & -1 &2 & -1 \\
0 & 0 & -1 & 2 \\
\end{bmatrix} ,\; \mathbf{u} =\begin{bmatrix}
u_1 \\ u_2 \\ u_3 \\ u_4
\end{bmatrix} $$
if we set $u_0 = u_5 = 0 $,  then $\frac{1}{\Delta^2}K\mathbf{u} = \mathbf{b}$ is the appropriate finite system. However it seems if we want arbitrary values we could solve instead
$\frac{1}{\Delta^2}K\mathbf{u} = \Big(\mathbf{b} + \frac{u_0}{\Delta^2}\mathbf{e}_1 + \frac{u_{n+1}}{\Delta^2}\mathbf{e}_n\Big)$. Is this the proper approach to the fix-fixed (Dirichlet) boundary conditions? What is the proper procedure for having free or Neumann/1st derivative boundary conditions instead? Practically, how does one construct the matrices to numerically solve a general 2nd order difference equation with Dirichlet or Neumann boundary conditions? How does this generalize if we embed a 2 dimensional problem into the matrix? 
 A: With Dirichlet conditions you can either build it into the right hand side of the equation, or you can introduce a trivial equation for those conditions as actual variables. It will be the same.
With Neumann or Neumann-like (e.g. Robin) conditions you want to be careful that you discretize the boundary derivative in a way which has the same order as your discretization of the interior derivatives. One way to do this with finite differences is to use "ghost points". I confess that this is rather hard to motivate within the finite difference framework but it gives results that are much like those you get in the finite element framework. 
The point here is that the values $u_0,u_{n+1}$ are thought of as being variables, and there is an equation involving the "interior" derivative at these points. You then add additional points $u_{-1},u_{n+2}$ and choose their values to satisfy the derivative condition. For example, with homogeneous Neumann conditions you take the approximation of the derivative at the first point to be $\frac{u_1-u_{-1}}{2\Delta}$ and set that to zero so that $u_{-1}=u_1$. Same for $u_{n+2}$. This then actually enters into the problem because you have an equation for the second derivative at the leftmost point: $u_{-1}-2u_0+u_1=b_0$, and $u_{-1}=u_1$ so this is $-2u_0+2u_1=b_0$.
