Updated sequential (non-looping) solver to include gamma*epsilon term.
Note that jmp = [u'] is now apparently not zero.