  We study a family of double confluent Heun equations of the form $\mathcal
E=0$, where $\mathcal L=\mathcal L_{\lambda,\mu,n}$ is a family of differential
operators of order two. They depend on complex parameters $\lambda$, $\mu$,
$n$. Its restriction to real parameter domain $\lambda+\mu^2>0$ is a
linearization of the family of nonlinear equations on two-torus modeling the
Josephson effect in superconductivity. We show that for every $b,n\in\mathbb C$
satisfying a certain "non-resonance condition" and every $\lambda,\mu\in\mathbb
C$, $\mu\neq0$ there exists an entire function $f_{\pm}:\mathbb C\to\mathbb C$
(unique up to constant factor) such that $z^{-b}\mathcal L(z^b
f_{\pm}(z^{\pm1}))=d_{0\pm}+d_{1\pm}z$ for some $d_{0\pm},d_{1\pm}\in\mathbb
C$. The constants $d_{j,\pm}$ are expressed as functions of the parameters.
This result has several applications. First of all, it gives the description of
those $\lambda$, $\mu$, $n$, $b$ for which the monodromy of the Heun equation
has eigenvalue $e^{2\pi i b}$. It also describes those $\lambda$, $\mu$, $n$
for which the monodromy is parabolic: has multiple eigenvalue. We consider the
rotation number $\rho$ of the dynamical system on two-torus as a function of
parameters restricted to a surface $\lambda+\mu^2=const$. The phase-lock areas
are its level sets having non-empty interiors. For general families of
dynamical systems the problem to describe the boundaries of the phase-lock
areas is known to be very complicated. Here we include the results in this
direction obtained by methods of complex variables. In our case the phase-lock
areas exist only for integer rotation numbers (quantization effect). The result
on parabolic monodromy implies the description of the union of their boundaries
by an explicit functional equation. For every $\theta\notin\mathbb Z$ we get a
description of the set $\{\rho\equiv\pm\theta(mod2\mathbb Z)\}$.
