The routine would not put elements at the right place in the residuals
vector (it would use variable indices instead of equation indices).

Also change the routine so that it computes a residual on evaluated
equations (instead of systematically returning zero for those).

Finally, simplify resid.m by calling this routine instead of embedding a
similar code.

Ref. #1823