Thanks for your points. I thought it would be cleaner if I pulled the below into a subroutine?
vec = (0.0d0, 0.0d0)
vec(1) = CMPLX(0.0d0,step*beta)
vec(numEq - 1) = (1.0d0, 0.0d0)
mat = (0.0d0,0.0d0)
mat(1,1) = (-1.0d0,0.0d0)
mat(1,2) = (1.0d0,0.0d0)
mat(1,numEq-1) =...