not really pretty but has all required outputs

The issue with GSS was due to the fact that the output is a Matrix with Boundary Values included, while my direct solvers deliver vectors without BV (that are later transformed in the PlotResults function). Everything works now. This version has Saumi's GSS, but Adam's works fine as well.