Let O subset of R-d be a bounded domain of class C-1,C-1. In L-2(O; C-n), we consider a selfadjoint second-order matrix elliptic differential operator B-N,B-epsilon,B- 0 < epsilon <= 1, under the Neumann boundary condition. The principal part of this operator is given in a factorized form. The operator includes first-order and zero-order terms. The coefficients of the operator B-N,B-epsilon are periodic and depend on x/epsilon. We study the generalized resolvent (B-N,B-epsilon - zeta Q(0)(center dot/epsilon))(-1), where Q(0) is a periodic bounded and positive definite matrix-valued function, and zeta is a complex parameter. We obtain approximations of the generalized resolvent in the operator norm in L-2(O;C-n) and in the norm of operators acting from L-2(O;C-n) to the Sobolev class H-1(O;C-n), with two-parametric (with respect to epsilon and zeta) error estimates. The results are applied to study the behavior of solutions of the initial boundary value problem with the Neumann condition for the parabolic equation Q(0)(x/epsilon)partial derivative(t)u(epsilon)(x,t)=-(B(N,epsilon)u(epsilon))(x,t) in the cylinder Ox(0,T), where 0< T <= infinity.