  Let $\mathcal{O}\subset\mathbb{R}^d$ be a bounded domain of class $C^{2p}$.
In $L_2(\mathcal{O};\mathbb{C}^n)$, we study a selfadjoint strongly elliptic
operator $A_{N,\varepsilon}$ of order $2p$ given by the expression $b({\mathbf
D})^* g({\mathbf x}/\varepsilon) b({\mathbf D})$, $\varepsilon >0$, with the
Neumann boundary conditions. Here $g({\mathbf x})$ is a bounded and positive
definite $(m\times m)$-matrix-valued function in ${\mathbb R}^d$, periodic with
respect to some lattice; $b({\mathbf D})=\sum_{|\alpha|=p} b_\alpha {\mathbf
D}^\alpha$ is a differential operator of order $p$ with constant coefficients;
$b_\alpha$ are constant $(m\times n)$-matrices. It is assumed that $m\geqslant
n$ and that the symbol $b({\boldsymbol \xi})$ has maximal rank for any $0 \ne
{\boldsymbol \xi}\in {\mathbb C}^d$. We find approximations for the resolvent
$\left(A_{N,\varepsilon}-\zeta I \right)^{-1}$ in the
$L_2(\mathcal{O};\mathbb{C}^n)$-operator norm and in the norm of operators
acting from $L_2(\mathcal{O};\mathbb{C}^n)$ to the Sobolev space
$H^p(\mathcal{O};\mathbb{C}^n)$, with error estimates depending on
$\varepsilon$ and $\zeta$.
