Programming interface for boundary conditions applied with blocks

[jlperla]

Point out in docs we do not implement the extension operator

• [x] pass in a list of operators with boundary already applied (come in as banded matrices or tridiagonal) and a Q matrix. It spits out the combined operator. I suggest the following interface

  bc = (Reflecting(), Reflecting())
  L1_bc = I * ρ - (μ1 * L₁₋bc(x̄, bc) - σ^2 / 2 * L₂bc(x̄, bc))
  L2_bc = I * ρ - (μ2 * L₁₋bc(x̄, bc) - σ^2 / 2 * L₂bc(x̄, bc))
  Q = [-0.05  0.05
        0.01 -0.01]
  L = jointoperator_bc( (L1_bc, L2_bc), Q)

  in there jointoperator_bc would check that the length of the tuple of passed in operators matches the size of the square Q matrix. (Resolved by 4df5fb1)

• [x] Make sure that you are using the best structure possible. Is https://github.com/JuliaMatrices/BlockBandedMatrices.jl the best? This isn't so trivial if the users pass in multiple types for the setups (e.g. a diffusion and a jump diffusion separately, where the first is trdiagonal and the second is banded or sparse). When you test, you may need to convert to banded matrices prior to stackign. (Resolved by 4df5fb1)
• [x] Add in unit tests for a 1x1 and 2x2 examples where you carefully check the matrices that pop out (Resolved by 4df5fb1)
• [x] Do a bigger 3x3 example with asset pricing and just check some numbers as a regression test (e.g. the asset pricing example we did before). (Resolved by 76393ed)
• [x] Test with a diffusion for the 2 states as well as a jump diffusion for the 2 states. (Resolved by 4df5fb1)
• [x] Make sure you can turn the block banded into sparse. (Run sparse(BANDED_BLOCK_MATRIX))

In all cases, be careful to verify that the types are what you are expecting. You may want to even do it with a check on returned type in the test.