BVP:
-u''(x) = f(x) = 2x + 1
u(0)    = g_0  = 3
u'(1)   = g_1  = -1/2
----------------------------
output:

creating equidistant mesh of (0, 1) with 20 elements...
stiffness matrix = 
TridiagonalMatrix(diag=[20, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 20], lower=[0, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20], upper=[0, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20])
load vector = 
[60, 60.055, 0.06, 0.065, 0.07, 0.075, 0.08, 0.085, 0.09, 0.095, 0.1, 0.105, 0.11, 0.115, 0.12, 0.125, 0.13, 0.135, 0.14, 0.145, -0.425]
system solved
u = [3, 3.07375, 3.14475, 3.21275, 3.2775, 3.33875, 3.39625, 3.44975, 3.499, 3.54375, 3.58375, 3.61875, 3.6485, 3.67275, 3.69125, 3.70375, 3.71, 3.70975, 3.70275, 3.68875, 3.6675]