\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 153, pp. 1--21.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2017/153\hfil Product-type system of difference equations]
{Two-dimensional product-type systems of difference equations of
 delay-type (2,2,1,2)}

\author[S. Stevi\'c \hfil EJDE-2017/153\hfilneg]
{Stevo Stevi\'c}

\address{Stevo Stevi\'c \newline
 Mathematical Institute of the Serbian Academy of Sciences,
Knez Mihailova 36/III,
11000 Beograd, Serbia. \newline
Operator Theory and  Applications Research Group,
Department of Mathematics,
King Abdulaziz University,
P.O. Box 80203, Jeddah 21589, Saudi Arabia}
\email{sstevic@ptt.rs}

\dedicatory{Communicated by Vicentiu Radulescu}

\thanks{Submitted April 21, 2017. Published June 27, 2017.}
\subjclass[2010]{39A20, 39A45}
\keywords{System of difference equations;  product-type system; 
\hfill\break\indent solvable in closed form}

\begin{abstract}
 We prove that the following class of systems of difference equations
 is solvable in closed form:
 $$
 z_{n+1}=\alpha z_{n-1}^aw_n^b,\quad
 w_{n+1}=\beta w_{n-1}^cz_{n-1}^d,\quad n\in\mathbb{N}_0,
 $$
 where
 $a, b, c, d\in\mathbb{Z}$, $\alpha, \beta, z_{-1}, z_0, w_{-1},
  w_0\in\mathbb{C}\setminus\{0\}$.
 We present formulas for its solutions in all the cases. The most complex
 formulas are presented in terms of the zeros of three different associated
 polynomials to the systems corresponding to the cases $a=0$, $c=0$ and
 $abcd\ne 0$, respectively, which on the other hand depend on some of
 parameters $a, b, c, d$.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\section{Introduction}

There has been a considerable interest in difference equations and systems
of difference equations (see, for example, 
\cite{ra,am1,amc218-sde, fp1,i5,j,kk,k,ll,jm,mk,p3-ejqtde,ps1,ps2,ps3,ps4,ps-jdea2011,
psh,pst2,rv,mr,sps,ps-cana,ejde-fodce,ejqtde-newm,ana2,ana1,
ejqtde-4th-2017,ejde-2015,jia2015-pts,ejde-2016,ade-16-ptdce,ejqtde-st2016}).
Books \cite{ra,j,k,ll,mk}, 
among others, present some classical methods for solving some classes of 
the equations and systems. The topic has re-attracted some attention
 in the last decade (see, for example,
 \cite{am1,amc218-sde,pst2,ejde-fodce, ejqtde-newm,ana2,ana1,ejqtde-4th-2017,
ejde-2015,jia2015-pts,ejde-2016,ade-16-ptdce,ejqtde-st2016}
and numerous related references therein). Many nonlinear equations and 
systems which have appeared in the literature recently were solved 
by transforming them to the classical solvable ones, by using some 
suitable changes of variables 
(see, for example, \cite{amc218-sde,pst2,ana1}, as well as the 
related references therein). On the other hand, almost two decades ago
 Papaschinopoulos and Schinas have started investigating some concrete
 systems of difference equations (\cite{ps1,ps2,ps3}),
 which motivated other experts to investigate some related ones
 (\cite{amc218-sde,fp1,p3-ejqtde,ps4,ps-jdea2011,psh,sps,ps-cana,ejde-fodce,
ejqtde-newm,ana2,ejqtde-4th-2017,ejde-2015,jia2015-pts,ejde-2016,ejde-2016,
ade-16-ptdce,ejqtde-st2016}).
 Majority of the systems investigated therein belong to the class of 
symmetric or close to symmetric systems. Namely, special cases of the 
following symmetric two-dimensional system of difference equations
$$
x_n=f(x_{n-k}, y_{n-l}),\quad y_n=f(y_{n-k}, x_{n-l}),
$$
$n\in\mathbb{N}_0$, have been studied a lot.

Many recent papers on difference equations and systems study only their 
positive solutions. One of the reasons for this is the fact that numerous
 difference equations and systems present some population models. 
Some systems of difference equations include product-type ones, which are
 obtained for some special values of its coefficients and/or parameters. 
It is well-known how product-type systems with positive initial values and 
coefficients are solved. Studying non-positive solutions is a more 
interesting problem. We started investigating complex-valued product-type 
systems in \cite{ejde-2015}. Another system was investigated
in \cite{jia2015-pts}, where some further basic steps concerning solvability 
of such systems were presented. Having published \cite{ana1} we came up with
 an idea of studying product-type systems which, unlike the ones in
 \cite{ejde-2015} and \cite{jia2015-pts}, have some multipliers.
The system in \cite{ejde-fodce} was the first complex-valued product-type 
system with multipliers that we have studied. 
Papers \cite{ejde-fodce,ejde-2015,jia2015-pts}, 
suggested us to study the solvability of the following two-dimensional 
product-type system
\begin{equation}
 z_n=\alpha z_{n-k}^aw_{n-l}^b,\quad w_n=\beta w_{n-m}^cz_{n-s}^d,\quad
n\in\mathbb{N},\label{gs}
\end{equation}
where $a, b, c, d, k, l, m, s\in\mathbb{N}$ and
$\alpha, \beta\in\mathbb{C}$. We call system \eqref{gs} two-dimensional
product-type system of {\it delay-type} $(k, m, l, s)$.

The main problem, which is a relatively big project, is to present 
closed-form formulas for the solvable systems of the type (this was 
not done for the systems in \cite{jia2015-pts} and \cite{ade-16-ptdce}). 
The systems in \cite{ana2} and \cite{ejqtde-st2016} were more complex 
than the one in \cite{ejde-fodce}, and non-trivial analyses of the form 
of their solutions in terms of the initial values and especially parameters 
were needed. The structures of the solutions for the system in \cite{ejde-2016}
 is simpler, so the corresponding analysis was also simpler. 
Another approach in dealing with the solvability problem for product-type 
systems can be found in \cite{ejqtde-newm}. To each system a few polynomials 
can be associated for determining the form of its solutions. 
A fourth order one helped in presenting formulas for solutions to the 
system in \cite{ejqtde-4th-2017}.


The solvability of system \eqref{gs} of delay-type $(2,2,1,2)$, that is, of 
\begin{equation}
z_{n+1}=\alpha z_{n-1}^aw_n^b,\quad
w_{n+1}=\beta w_{n-1}^cz_{n-1}^d,\quad n\in\mathbb{N}_0,\label{ms1}
\end{equation}
where $a,b,c,d\in\mathbb{Z}$, $\alpha,\beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}$,
is investigated in this paper, continuing the project in
 \cite{ejde-fodce,ejqtde-newm,ana2,ejqtde-4th-2017,ejde-2015,jia2015-pts,
ejde-2016,ade-16-ptdce,ejqtde-st2016}.

Note that if some of the quantities $\alpha, \beta, w_{-1}, w_0, z_{-1}, z_0$ 
are zero, then the solutions are either not defined or trivial, so not 
interesting. Hence, the case is excluded from further considerations. 
If $m,n\in\mathbb{Z}$, then $\overline{m,n}=\{j\in\mathbb{Z}$, $m\le j\le n\}$, 
while $\sum_{j=n}^{n-1}f_j=0$ for each $n\in\mathbb{Z}$.

\section{Auxiliary results}

The following set of lemmas presents standard tools for dealing with 
the problem of solvability of product-type systems with small delays 
$k, l, m, s$. For the first lemma, which is  classical one, see, 
for example, \cite{f, k, jia2015-pts}.

\begin{lemma} \label{lem1} 
 Let
$$
p_k(t)=c_k\prod_{j=1}^k(t-t_j),
$$
$c_k\ne 0$ and $t_i\ne t_j$, $i\ne j$. Then
\begin{gather*}
\sum_{j=1}^k\frac{t_j^s}{p'_k(t_j)}=0,\quad 0\le s\le k-2, \\
\sum_{j=1}^k\frac{t_j^{k-1}}{p'_k(t_j)}=\frac1{c_k}.
\end{gather*}
\end{lemma}

For the second lemma see, for example, \cite{k, mk} (or \cite{ana2} 
for a general method for calculating the sums in the lemma).

\begin{lemma} \label{lem2}  Let
$$
s_n^{(i)}(z)=\sum_{j=1}^nj^iz^{j-1},\quad n\in\mathbb{N},
$$
where $i\in\mathbb{N}_0$ and $z\in\mathbb{C}$.
Then
\begin{gather*}
s_n^{(0)}(z)=\frac{1-z^n}{1-z},\\
s_n^{(1)}(z)=\frac{1-(n+1)z^n+nz^{n+1}}{(1-z)^2},\\
s_n^{(2)}(z)=\frac{1+z-(n+1)^2z^n+(2n^2+2n-1)z^{n+1}-n^2z^{n+2}}{(1-z)^3},\\
s_n^{(3)}(z)=\frac{n^3z^n(z-1)^3-3n^2z^n(z-1)^2+3nz^n(z^2-1)
-(z^n-1)(z^2+4z+1)}{(1-z)^4},
\end{gather*}
for every $z\in\mathbb{C}\setminus\{1\}$ and $n\in\mathbb{N}$.
\end{lemma}

The results in \cite{r} are suitably reformulated in the following lemma.

\begin{lemma} \label{lem3}  Let
\begin{gather*}
P_4(t)=t^4+bt^3+ct^2+dt+e,\quad \Delta_0=c^2-3bd+12e, \\
 \Delta_1=2c^3-9bcd+27b^2e+27d^2-72ce,\quad
\Delta=\frac1{27}(4\Delta_0^3-\Delta_1^2), \\
P=8c-3b^2,\quad Q=b^3+8d-4bc,\quad D=64e-16c^2+16b^2c-16bd-3b^4.
\end{gather*}
\begin{itemize}
\item[(a)] If $\Delta<0$, then two zeros of $P_4$ are real and different, 
 and two are complex conjugate;
\item[(b)] If $\Delta>0$, then all the zeros of $P_4$ are real or none is.
  More precisely,
    \begin{itemize}
    \item[1] if $P<0$ and $D<0$, then all four zeros of $P_4$ are real and different;
    \item[2] if $P>0$ or $D>0$, then there are two pairs of complex conjugate zeros of $P_4$.
    \end{itemize}
\item[(c)] If $\Delta=0$, then and only then $P_4$ has a multiple zero. The following cases can occur:
    \begin{itemize}
    \item[1] if $P<0$, $D<0$ and $\Delta_0\ne 0$, then two zeros of $P_4$ are
  real and equal and two are real and simple;
    \item[2] if $D>0$ or ($P>0$ and ($D\ne 0$ or $Q\ne 0$)), then two zeros 
 of $P_4$ are real and equal and two are complex conjugate;
    \item[3] if $\Delta_0=0$ and $D\ne 0$, there is a triple zero of $P_4$ and 
 one simple, all real;
    \item[4] if $D=0$, then
        \begin{itemize}
        \item[4.1] if $P<0$ there are two double real zeros of $P_4$;
        \item[4.2] if $P>0$ and $Q=0$ there are two double complex conjugate 
 zeros of $P_4$;
        \item[4.3] if $\Delta_0=0$, then all four zeros of $P_4$ are real and 
 equal to $-b/4$.
        \end{itemize}
    \end{itemize}
\end{itemize}
\end{lemma}


\section{Main results}


This section contains our main results. Before formulating them we 
give a list of a few members of sequences $z_n$ and $w_n$ that we use 
in proofs of the results
\begin{equation}
\begin{gathered}
z_1=\alpha z_{-1}^aw_0^b,\quad z_2=\alpha\beta^bw_{-1}^{bc}z_{-1}^{bd}z_0^a\\
w_1=\beta w_{-1}^c z_{-1}^d,\quad w_2=\beta w_0^c z_0^d.
\end{gathered}\label{z}
\end{equation}

\begin{theorem} \label{thm1} 
Assume that $a,c,d\in\mathbb{Z}$, $b=0$, $\alpha,\beta, z_{-1}, z_0, w_{-1},
w_0\in\mathbb{C}\setminus\{0\}$. Then system
\eqref{ms1} is solvable in closed form.
\end{theorem}

\begin{proof} 
From the condition $b=0$, we have
\begin{equation}
z_{n+1}=\alpha z_{n-1}^a,\quad
w_{n+1}=\beta w_{n-1}^cz_{n-1}^d,\quad n\in\mathbb{N}_0.\label{ms2}
\end{equation}
The first equation in \eqref{ms2} yields
\begin{equation}
 z_{2n+i}=\alpha^{\sum_{j=0}^{n-1}a^j}z_i^{a^n},\label{c1}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$, from which it follows that
\begin{equation}
 z_{2n+i}=\alpha^{\frac{1-a^n}{1-a}}z_i^{a^n},\label{c2a}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$, when $a\ne 1$, and
\begin{equation}
 z_{2n+i}=\alpha^nz_i,\label{c3}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$, when $a=1$.

Using \eqref{c1} in the second equation in \eqref{ms2} we obtain
\begin{equation}
 w_{2n+i}=\beta \alpha^{d\sum_{j=0}^{n-2}a^j}z_i^{da^{n-1}}w_{2(n-1)+i}^c,\label{c5}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.

From \eqref{c5} and by induction it is proved that
\begin{equation}
 w_{2n+i}=\beta^{\sum_{l=0}^{n-1}c^l}\alpha^{d\sum_{l=0}^{n-1}c^l\sum_{j=0}^{n-l-2}a^j}
z_i^{d\sum_{l=0}^{n-1}c^la^{n-l-1}}w_i^{c^n},\label{c6a}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.
\smallskip

\noindent\textbf{Case $a\ne 1\ne c\ne a$.}
 From \eqref{c6a} and by Lemma \ref{lem2}, in this case, we have
\begin{equation}
\begin{aligned}
w_{2n+i}=&\beta^{\sum_{l=0}^{n-1}c^l}\alpha^{d\sum_{l=0}^{n-1}c^l
 \sum_{j=0}^{n-l-2}a^j}z_i^{d\sum_{l=0}^{n-1}c^la^{n-l-1}}w_i^{c^n}\\
=&\beta^{\frac{1-c^n}{1-c}}\alpha^{d\sum_{l=0}^{n-1}c^l
 \frac{1-a^{n-l-1}}{1-a}}z_i^{d\frac{a^n-c^n}{a-c}}w_i^{c^n}\\
=&\beta^{\frac{1-c^n}{1-c}}\alpha^{\frac{d}{1-a}(\frac{1-c^n}{1-c}
 -\frac{a^n-c^n}{a-c})}z_i^{d\frac{a^n-c^n}{a-c}}w_i^{c^n}\\
=&\beta^{\frac{1-c^n}{1-c}}\alpha^{\frac{d(a-c+(1-a)c^n
 +(c-1)a^n)}{(1-a)(1-c)(a-c)}}z_i^{d\frac{a^n-c^n}{a-c}}w_i^{c^n},
\end{aligned}\label{c6}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.
\smallskip


\noindent\textbf{Case $a=c\ne 1$.} 
From \eqref{c6a} and by Lemma \ref{lem2}, in this case, we have
\begin{equation}
\begin{aligned}
w_{2n+i}=&\beta^{\sum_{l=0}^{n-1}a^l}\alpha^{d\sum_{l=0}^{n-1}a^l
 \sum_{j=0}^{n-l-2}a^j}z_i^{d\sum_{l=0}^{n-1}a^la^{n-l-1}}w_i^{a^n}\\
=&\beta^{\frac{1-a^n}{1-a}}\alpha^{d\sum_{l=0}^{n-1}a^l
 \frac{1-a^{n-l-1}}{1-a}}z_i^{dna^{n-1}}w_i^{a^n}\\
=&\beta^{\frac{1-a^n}{1-a}}\alpha^{\frac{d}{1-a}(\frac{1-a^n}{1-a}
 -na^{n-1})}z_i^{dna^{n-1}}w_i^{a^n}\\
=&\beta^{\frac{1-a^n}{1-a}}\alpha^{\frac{d((n-1)a^n-na^{n-1}+1)}
 {(1-a)^2}}z_i^{dna^{n-1}}w_i^{a^n},
\end{aligned} \label{c7}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.
\smallskip

\noindent\textbf{Case $a\ne 1=c$.} 
From \eqref{c6a} and by Lemma \ref{lem2}, in this case, we have
\begin{equation}
\begin{aligned}
w_{2n+i}=&\beta^{\sum_{l=0}^{n-1}1}\alpha^{d\sum_{l=0}^{n-1}
 \sum_{j=0}^{n-l-2}a^j}z_i^{d\sum_{l=0}^{n-1}a^{n-l-1}}w_i\\
=&\beta^n\alpha^{d\sum_{l=0}^{n-1}\frac{1-a^{n-l-1}}{1-a}}z_i^{d\frac{a^n-1}{a-1}}w_i\\
=&\beta^n\alpha^{\frac{d}{1-a}(n-\frac{a^n-1}{a-1})}z_i^{d\frac{a^n-1}{a-1}}w_i\\
=&\beta^n\alpha^{d\frac{a^n-an+n-1}{(1-a)^2}}z_i^{d\frac{a^n-1}{a-1}}w_i,
\end{aligned}\label{c8}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0.$
\smallskip

\noindent\textbf{Case $a=1\ne c$.} 
From \eqref{c6a} and by Lemma \ref{lem2}, in this case, we have
\begin{equation}
\begin{aligned}
w_{2n+i}=&\beta^{\sum_{l=0}^{n-1}c^l}\alpha^{d\sum_{l=0}^{n-1}c^l(n-l-1)}z_i
 ^{d\sum_{l=0}^{n-1}c^l}w_i^{c^n}\\
=&\beta^{\frac{1-c^n}{1-c}}\alpha^{d((n-1)\frac{1-c^n}{1-c}
 -c\frac{1-nc^{n-1}+(n-1)c^n}{(1-c)^2}}z_i^{d\frac{1-c^n}{1-c}}w_i^{c^n}\\
=&\beta^{\frac{1-c^n}{1-c}}\alpha^{d\frac{c^n-nc+n-1}{(1-c)^2}}z_i
 ^{d\frac{1-c^n}{1-c}}w_i^{c^n},
\end{aligned} \label{c9}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.
\smallskip

\noindent\textbf{Case $a=c=1$.} 
From \eqref{c6a} and by a well-known formula, in this case, we have
\begin{equation}
\begin{aligned}
w_{2n+i}=&\beta^{\sum_{l=0}^{n-1}1}\alpha^{d\sum_{l=0}^{n-1}(n-l-1)}z_i
 ^{d\sum_{l=0}^{n-1}1}w_i\\
=&\beta^n\alpha^{d\frac{(n-1)n}2}z_i^{dn}w_i,
\end{aligned} \label{c10}
\end{equation}
for $n\in\mathbb{N}$ and $i=-1,0$.

From all above mentioned the theorem follows.
\end{proof}

\begin{corollary} \label{coro1} 
 Assume that $a,c,d\in\mathbb{Z}$, $b=0$, $\alpha, \beta, z_{-1}, z_0, 
w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. Then the following statements are true.
\begin{itemize}
\item[(a)]  If $a\ne 1\ne c\ne a$, then the general solution to system
\eqref{ms1} is given by formulas \eqref{c2a} and \eqref{c6}.

\item[(b)] If $a=c\ne 1$, then the general solution to system
\eqref{ms1} is given by formulas \eqref{c2a} and \eqref{c7}.

\item[(c)] If $a\ne 1=c$, then the general solution to system
\eqref{ms1} is given by formulas \eqref{c2a} and \eqref{c8}.

\item[(d)] If $a=1\ne c$, then the general solution to system
\eqref{ms1} is given by formulas \eqref{c3} and \eqref{c9}.

\item[(e)] If $a=c=1$, then the general solution to system
\eqref{ms1} is given by formulas \eqref{c3} and \eqref{c10}.

\end{itemize}
\end{corollary}


\begin{theorem} \label{thm2}
 Assume that $a,b,c\in\mathbb{Z}$, $d=0$, 
$\alpha,\beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. Then system
\eqref{ms1} is solvable in closed form.
\end{theorem}

\begin{proof} 
Using the condition $d=0$ in \eqref{ms1} we have
\begin{equation}
z_{n+1}=\alpha z_{n-1}^a w_n^b,\quad
w_{n+1}=\beta w_{n-1}^c,\quad n\in\mathbb{N}_0.\label{ms3}
\end{equation}
This system  corresponds  the one in \cite[Theorem 2.3]{ejde-2016},
 where parameters $a$ and $b$, as well as $c$ and $d$ are interchanged,
respectively. Hence, all the formulas for solutions to \eqref{ms3}
are  obtained from \cite[Theorem 2.2 and Corollary 2.2]{ejde-2016}.
\end{proof}


\begin{theorem} \label{thm3}
 Assume that $b, c, d\in\mathbb{Z}$, $bd\ne 0$, $a=0$, 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. Then system
\eqref{ms1} is solvable in closed form.
\end{theorem} 

\begin{proof}
 Using the condition $a=0$ in \eqref{ms1}, we have
\begin{equation}
 z_{n+1}=\alpha w_n^b,\quad
w_{n+1}=\beta w_{n-1}^cz_{n-1}^d,\quad n\in\mathbb{N}_0.\label{q1}
\end{equation}
From \eqref{q1} we easily get
\begin{equation}
 w_{n+1}=\alpha^d\beta w_{n-1}^cw_{n-2}^{bd},\quad n\ge 2.\label{q5w}
\end{equation}
Let $\mu=\alpha^d\beta$,
\begin{equation}
 a_1=0,\quad b_1=c,\quad c_1=bd,\quad y_1=1.\label{icq}
\end{equation}
Then
\begin{equation}
w_{n+1}=\mu^{y_1}w_n^{a_1}w_{n-1}^{b_1}w_{n-2}^{c_1},\quad
n\ge 2.\label{d6q}
\end{equation}
Similarly, \eqref{d6q} implies
\begin{equation}
\begin{aligned}
w_{n+1}
&=\mu^{y_1}(\mu w_{n-1}^{a_1}w_{n-2}^{b_1}w_{n-3}^{c_1})^{a_1}w_{n-1}^{b_1}w_{n-2}^{c_1},\\
&=\mu^{y_1+a_1}w_{n-1}^{a_1a_1+b_1}w_{n-2}^{b_1a_1+c_1}w_{n-3}^{c_1a_1}\\
&=\mu^{y_2}w_{n-1}^{a_2}w_{n-2}^{b_2}w_{n-3}^{c_2},
\end{aligned} \label{d8q}
\end{equation}
for $n\ge 3$, where
\begin{equation}
a_2:=a_1a_1+b_1,\quad b_2:=b_1a_1+c_1,\quad c_2:=c_1a_1,\quad
y_2:=y_1+a_1.\label{d85q}
\end{equation}
Assume
\begin{equation}
w_{n+1}=\mu^{y_k} w_{n+1-k}^{a_k}w_{n-k}^{b_k}w_{n-k-1}^{c_k},\label{d9q}
\end{equation}
for a $k\ge 2$ and all $n\ge k+1$, and
\begin{gather}
a_k=a_1a_{k-1}+b_{k-1},\quad b_k=b_1a_{k-1}+c_{k-1},\quad c_k=c_1a_{k-1},\label{d10q}\\
y_k=y_{k-1}+a_{k-1}.\label{d10aq}
\end{gather}


If we replace $n$ by $n-k$ in \eqref{d6q}, and employ it in \eqref{d9q}, we obtain
\begin{equation}
\begin{aligned}
w_{n+1}
&=\mu^{y_k}(\mu w_{n-k}^{a_1}w_{n-k-1}^{b_1}w_{n-k-2}^{c_1})^{a_k}w_{n-k}^{b_k}w_{n-k-1}^{c_k}\\
&=\mu^{y_k+a_k}w_{n-k}^{a_1a_k+b_k}w_{n-k-1}^{b_1a_k+c_k}w_{n-k-2}^{c_1a_k}\\
&=\mu^{y_{k+1}}w_{n-k}^{a_{k+1}}w_{n-k-1}^{b_{k+1}}w_{n-k-2}^{c_{k+1}},
\end{aligned} \label{d11q}
\end{equation}
for $n\ge k+2$, where
\begin{equation}
a_{k+1}:=a_1a_k+b_k,\quad b_{k+1}:=b_1a_k+c_k,\quad
c_{k+1}:=c_1a_k,\quad y_{k+1}:=y_k+a_k.\label{d12q}
\end{equation}
Equalities \eqref{d8q}, \eqref{d85q}, \eqref{d11q}, \eqref{d12q} along with the
induction show that \eqref{d9q}-\eqref{d10aq} hold.

Setting $k=n-1$ in \eqref{d9q} and using \eqref{z}, \eqref{d10q} and \eqref{d10aq},
 we have
\begin{equation}
\begin{aligned}
w_{n+1}&=\mu^{y_{n-1}}w_2^{a_{n-1}}w_1^{b_{n-1}}w_0^{c_{n-1}}\\
&=(\alpha^d\beta)^{y_{n-1}}(\beta w_0^c z_0^d)^{a_{n-1}}(\beta w_{-1}^c
  z_{-1}^d)^{b_{n-1}}w_0^{c_{n-1}}\\
&=\alpha^{dy_{n-1}}\beta^{y_{n-1}+a_{n-1}}z_{-1}^{db_{n-1}}z_0^{da_{n-1}}
 w_{-1}^{cb_{n-1}}w_0^{ca_{n-1}+c_{n-1}}\\
&=\alpha^{dy_{n-1}}\beta^{y_n}z_{-1}^{da_n}z_0^{da_{n-1}}w_{-1}^{ca_n}
 w_0^{a_{n+1}},
\end{aligned} \label{d19q}
\end{equation}
for $n\geq 2$.


From \eqref{d10q} we obtain
 \begin{equation}
a_k=b_1a_{k-2}+c_1a_{k-3},\quad\text{for } k\ge 4.\label{d17q}
\end{equation}
Since $c_1=bd\ne 0$, from \eqref{d17q} it follows that
$$
a_{k-3}=\frac{a_k-b_1a_{k-2}}{c_1}.
$$
Using this equality, $a_j$-s can be calculated for $j\le 0$, and,
as in \cite{ana2} and \cite{ejqtde-st2016}, we obtain
\begin{equation}
 a_{-1}=a_{-2}=0,\quad a_0=1.\label{b1q}
\end{equation}
From \eqref{icq}, \eqref{d10aq} and \eqref{b1q}, it follows that
\begin{gather}
 y_{-2}=y_{-1}=y_0=0,\quad y_1=1,\label{ic3q} \\
 y_k=\sum_{j=0}^{k-1}a_j,\quad k\in\mathbb{N}.\label{ynq}
\end{gather}

Since the initial-value problem \eqref{d17q}-\eqref{b1q} is solvable, 
formula for $a_k$ can be found, from which along with \eqref{ynq} and 
by Lemma \ref{lem2} is obtained a formula for $y_k$ (the form of $a_k$ is well-known,
 see \cite{ra, mk}). Using the formulas for $a_k$ and $y_k$ in \eqref{d19q} 
we have a closed form formula for the solution to equation \eqref{q5w}.


Using \eqref{d19q} in the first equation in \eqref{q1}, we obtain
\begin{equation}
z_n
=\alpha^{1+bdy_{n-3}}\beta^{by_{n-2}}z_{-1}^{bda_{n-2}}z_0^{bda_{n-3}}w_{-1}^{bca_{n-2}}
w_0^{ba_{n-1}}, \label{f19q}
\end{equation}
for $n\in\mathbb{N}$.

It is not difficult to see that formulas \eqref{d19q} and \eqref{f19q} 
present a solution to system \eqref{q1}. This completes the proof.
\end{proof} 


Now we specify the form of the solutions to \eqref{ms1} when $a=0$ and $bd\ne 0$. 
The characteristic polynomial associated to \eqref{d17q} is
\begin{equation}
 p_3(\lambda)=\lambda^3-c\lambda-bd.\label{t27k}
\end{equation}
The zeros of \eqref{t27k} are
\begin{equation} \lambda_j=\frac{1}{3\sqrt[3]{2}}\Big(\varepsilon^j\sqrt[3]
{\Delta_1-\sqrt{\Delta_1^2-4\Delta_0^3}}
+\overline{\varepsilon}^j\sqrt[3]{\Delta_1+\sqrt{\Delta_1^2-4\Delta_0^3}}\Big),
\quad j=\overline{0,2},\label{s2g}
\end{equation}
where
\begin{equation}
\Delta_0=3c=:-3p\quad\text{and}\quad \Delta_1=27bd=:-27q,\label{de1g}
\end{equation}
and $\varepsilon^3=1$, $\varepsilon\ne1$ (\cite{f}).
\smallskip

\noindent\textbf{Zeros of $p_3$ are different and none of them is 1.}
In this case it must be $\Delta:=(\Delta_1^2-4\Delta_0^3)/27\ne 0$, that is, 
$27(bd)^2\ne 4c^3$.
 If also $c+bd\ne 1$, then 1 is not a zero of $p_3$. 
These conditions are satisfied, for example, if $c=bd\in\mathbb{N}$.
\smallskip

\noindent\textbf{Zeros of $p_3$ are different and one of them is 1.} 
Since $p_3(1)=0$, we have that $c+bd=1$.
Hence
$$
p_3(\lambda)=\lambda^3-c\lambda+c-1=(\lambda-1)(\lambda^2+\lambda+1-c),
$$
and the zeros of $p_3$ are
\begin{equation}
 \lambda_{1,2}=\frac{-1\pm\sqrt{4c-3}}2,\quad \lambda_3=1.\label{l5}
\end{equation}
Note also that it must be $c\ne 3$, from the condition $p_3'(1)=3-c\ne 0$.


Since in these two cases all the zeros of $p_3$ are different, 
the general solution to \eqref{d17q} must have the  form
\begin{equation}
a_n=\alpha_1\lambda_1^n+\alpha_2\lambda_2^n+\alpha_3\lambda_3^n,\quad
n\in\mathbb{N},\label{d51m}
\end{equation}
where $\alpha_i$, $i=\overline{1,3}$, are constants \cite{ra, mk}.

As it was explained, for example, in \cite{ana2}, the solution to 
\eqref{d17q} satisfying \eqref{b1q} is given by
\begin{equation}
\begin{aligned}
 a_n&=\sum_{j=1}^3 \frac{\lambda_j^{n+2}}{p_3'(\lambda_j)}\\
&=  \frac{\lambda_1^{n+2}}{(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)}
 + \frac{\lambda_2^{n+2}}{(\lambda_2-\lambda_1)(\lambda_2-\lambda_3)}
 + \frac{\lambda_3^{n+2}}{(\lambda_3-\lambda_1)(\lambda_3-\lambda_2)},
\end{aligned}\label{s1m}
\end{equation}
and the formula holds for every $n\in\mathbb{Z}$.

From \eqref{ynq} and \eqref{s1m}, it follows that
\begin{equation}
y_n=\sum_{k=0}^{n-1}\sum_{j=1}^3\frac{\lambda_j^{k+2}}{p_3'(\lambda_j)},
\quad n\in\mathbb{N}.\label{a27}
\end{equation}
The formula holds also for $n\ge -2$ (\cite{ana2}).

If $\lambda_j\ne 1$, $j=\overline{1,3}$, then from \eqref{a27} and Lemma \ref{lem2}, we obtain
\begin{equation}
\begin{aligned}
y_n&=\frac{\lambda_1^2(\lambda_1^n-1)}{(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)(\lambda_1-1)}
+\frac{\lambda_2^2(\lambda_2^n-1)}{(\lambda_2-\lambda_1)(\lambda_2-\lambda_3)(\lambda_2-1)}\\
&\quad +\frac{\lambda_3^2(\lambda_3^n-1)}{(\lambda_3-\lambda_1)(\lambda_3-\lambda_2)(\lambda_3-1)},
\end{aligned} \label{a28b}
\end{equation}
while if one of the zeros is 1, say $\lambda_3$, then $1\ne\lambda_1\ne \lambda_2\ne 1$,
 and we have
\begin{equation}
y_n=\frac{\lambda_1^2(\lambda_1^n-1)}{(\lambda_1-\lambda_2)(\lambda_1-1)^2}
+\frac{\lambda_2^2(\lambda_2^n-1)}{(\lambda_2-\lambda_1)(\lambda_2-1)^2}
+\frac{n}{(\lambda_1-1)(\lambda_2-1)}.\label{a29b}
\end{equation}
for $n\in\mathbb{N}$. Formulas \eqref{a28b} and \eqref{a29b} holds for every
$n\ge -2$ (\cite{ana2}).


\begin{corollary} \label{coro2} 
 Assume that $b,c,d\in\mathbb{Z}$, $a=0$, $bd\ne0$, 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$ and
 $\Delta\ne 0$. Then the following statements are true.
\begin{itemize}
\item[(a)] If $c+bd\ne 1$, then the general solution to
\eqref{ms1} is given by \eqref{d19q} and \eqref{f19q}, 
where $(a_n)_{n\ge-2}$ is given by
\eqref{s1m}, $(y_n)_{n\ge-2}$ is given by \eqref{a28b}, while $\lambda_j$-s, 
$j=\overline{1,3}$, are given by \eqref{s2g} and \eqref{de1g}.

\item[(b)] If $c+bd=1$ and $c\ne 3$, then $p_3$ has a unique zero equal to 1, 
say $\lambda_3$, and the general solution to
\eqref{ms1} is given by formulas \eqref{d19q} and \eqref{f19q}, where
 $(a_n)_{n\ge-2}$ is given by
\eqref{s1m} with $\lambda_3=1$, $(y_n)_{n\ge-2}$ is given by \eqref{a29b}, 
while $\lambda_j$-s, $j=\overline{1,3}$, are given by \eqref{l5}.

\end{itemize}
\end{corollary}
\smallskip

\noindent\textbf{One of the zeros is double.}
 In this case it must be $\Delta=0$, that is, $c^3=27(bd)^2/4.$ If $m$ is a 
double zero of $p_3$, then
$$
m^3-cm-bd=0\quad\text{and}\quad 3m^2-c=0.
$$
Hence
\begin{equation}
 p_3(\lambda)=\lambda^3-3m^2\lambda+2m^3=(\lambda-m)^2(\lambda+2m).\label{t28}
\end{equation}
Since $bd\ne 0$, we have $m\ne 0$. From this and \eqref{t28} we see that
$p_3$ cannot have a triple zero. It also cannot have a unique zero equal
to 1, otherwise we would have $2m=1$ and consequently $c=3/4\notin\mathbb{Z}$,
which is impossible. Note also that 1 is a double zero if and only if $c=3$
and $bd=-2$.

If $\lambda_1\ne\lambda_2=\lambda_3$, then the general solution to \eqref{d17q} has the  form
\begin{equation}
 a_n=\hat \alpha_1\lambda_1^n+(\hat \alpha_2+\hat \alpha_3n)\lambda_2^n,\quad
n\in\mathbb{N},\label{a15a}
\end{equation}
where $\hat\alpha_i$, $i=\overline{1,3}$, are constants.

The solution of the form in \eqref{a15a} satisfying \eqref{b1q} is
\begin{equation}
 a_n=\frac{\lambda_1^{n+2}+\big(\lambda_2-2\lambda_1+n(\lambda_2-\lambda_1)\big)\lambda_2^{n+1}}{(\lambda_2-\lambda_1)^2},
\label{a16}
\end{equation}
and also holds for every $n\ge -2$ (\cite{ana2}).

Since, in the case, $\Delta=0$, from \eqref{s2g}, we have
\begin{equation}
\lambda_1=\frac{2}{3\sqrt[3]{2}}\sqrt[3]{\Delta_1},\quad \lambda_{2,3}
=-\frac1{3\sqrt[3]{2}}\sqrt[3]{\Delta_1}.\label{w7}
\end{equation}
From \eqref{ynq} and \eqref{a16}, we have
\begin{equation}
y_n=\sum_{j=0}^{n-1}a_j=\sum_{j=0}^{n-1}\frac{\lambda_1^{j+2}
+\big(\lambda_2-2\lambda_1+j(\lambda_2-\lambda_1)\big)\lambda_2^{j+1}}{(\lambda_2-\lambda_1)^2},\quad
n\in\mathbb{N}.\label{a32}
\end{equation}
The formula holds also for $n\ge -2$ (\cite{ana2}).

From \eqref{a32} and by Lemma \ref{lem2}, it follows that
\begin{equation}
\begin{aligned}
y_n&=\frac{\lambda_1^2(\lambda_1^n-1)}{(\lambda_2-\lambda_1)^2(\lambda_1-1)}+
\frac{(\lambda_2-2\lambda_1)\lambda_2(\lambda_2^n-1)}{(\lambda_2-\lambda_1)^2(\lambda_2-1)}\\
&\quad +\frac{\lambda_2^2(1-n\lambda_2^{n-1}+(n-1)\lambda_2^n)}{(\lambda_2-\lambda_1)(\lambda_2-1)^2},\quad
 n\in\mathbb{N}.\label{a32f}
\end{aligned}
\end{equation}

If we assume that $\lambda_1\ne 1$ and $\lambda_2=\lambda_3=1$, then from \eqref{a32} 
it follows that
\begin{equation}
y_n=\frac{\lambda_1^2(\lambda_1^n-1)}{(\lambda_1-1)^3}+
\frac{(1-2\lambda_1)n}{(\lambda_1-1)^2}+\frac{(n-1)n}{2(1-\lambda_1)},\quad
n\in\mathbb{N}.\label{t29}
\end{equation}
Formulas \eqref{a32f} and \eqref{t29} hold also for $n\ge -2$ (\cite{ana2}).


\begin{corollary} \label{coro3}
 Assume that
$b,c,d\in\mathbb{Z}$, $a=0$, $bd\ne 0$,  $\alpha, \beta, z_{-1}, z_0, w_{-1},
w_0\in\mathbb{C}\setminus\{0\}$ and $\Delta=0$. Then the following statements
 are true.
\begin{itemize}
\item[(a)]  If $c+bd\ne 1$, then the general solution to
\eqref{ms1} is given by \eqref{d19q} and \eqref{f19q}, where 
$(a_n)_{n\ge-2}$ is given by
\eqref{a16}, $(y_n)_{n\ge-2}$ is given by \eqref{a32f}, while 
$\lambda_j$-s, $j=\overline{1,3}$, are given by \eqref{w7} and \eqref{de1g}.

\item[(b)] If $c=3$ and $bd=-2$, then two zeros of \eqref{t27k} are equal to 1,
 say, $\lambda_2$ and $\lambda_3$, and the general solution to system
\eqref{ms1} is given by \eqref{d19q} and \eqref{f19q},  where 
$(a_n)_{n\ge-2}$ is given by
\eqref{a16} with $\lambda_2=1$, $(y_n)_{n\ge-2}$ is given by \eqref{t29},
 while $\lambda_1=-2$.

\item[(c)] Polynomial \eqref{t27k} cannot have 1 as a simple zero.

\item[(d)] Polynomial \eqref{t27k} cannot have a triple zero.

\end{itemize}
\end{corollary}


\begin{theorem} \label{thm4} 
 Assume that $a,b,d\in\mathbb{Z}$, $bd\ne 0$, $c=0$, 
$\alpha, \beta, z_{-1}, z_0, w_0\in\mathbb{C}\setminus\{0\}$. Then system
\eqref{ms1} is solvable in closed form.
\end{theorem}

\begin{proof} 
The condition $c=0$ implies that
\begin{equation} z_{n+1}=\alpha z_{n-1}^aw_n^b,\quad
w_{n+1}=\beta z_{n-1}^d,\quad n\in\mathbb{N}_0.\label{r1a}
\end{equation}
From \eqref{r1a} we obtain
\begin{equation}
 z_{n+1}=\alpha\beta^b z_{n-1}^az_{n-2}^{bd},\quad n\in\mathbb{N}.\label{r5a}
\end{equation}
Let $\mu=\alpha\beta^b$,
\begin{equation}
a_1=0,\quad b_1=a,\quad c_1=bd,\quad y_1=1.\label{icz}
\end{equation}
Then
\begin{equation}
z_{n+1}=\mu^{y_1}z_n^{a_1}z_{n-1}^{b_1}z_{n-2}^{c_1},\quad
n\in\mathbb{N}.\label{d6z}
\end{equation}

As in the proof of Theorem \ref{thm3} we obtain
\begin{gather}
z_{n+1}=\mu^{y_k} z_{n+1-k}^{a_k}z_{n-k}^{b_k}z_{n-k-1}^{c_k},\label{d9z}\\
a_k=a_1a_{k-1}+b_{k-1},\quad b_k=b_1a_{k-1}+c_{k-1},
\quad c_k=c_1a_{k-1},\label{d10z}\\
y_k=y_{k-1}+a_{k-1},\label{d10az}
\end{gather}
for each $k\ge 2$ and all $n\ge k$.


Setting $k=n$ in \eqref{d9z} and using \eqref{z}, \eqref{d10z} and \eqref{d10az}, we have
\begin{align}
z_{n+1}&=\mu^{y_n}z_1^{a_n}z_0^{b_n}z_{-1}^{c_n}\\
&=(\alpha\beta^b)^{y_n}(\alpha z_{-1}^aw_0^b)^{a_n}z_0^{b_n}z_{-1}^{c_n}\\
&=\alpha^{y_n+a_n}\beta^{by_n}z_{-1}^{aa_n+c_n}z_0^{b_n}w_0^{ba_n}\\
&=\alpha^{y_{n+1}}\beta^{by_n}z_{-1}^{a_{n+2}}z_0^{a_{n+1}}w_0^{ba_n},\label{d19z}\end{align}
for $n\in\mathbb{N}.$

From \eqref{d10z},
 \begin{equation}
a_k=b_1a_{k-2}+c_1a_{k-3},\qquad\text{for } k\ge 4.\label{d17z}
\end{equation}
Since $c_1=bd\ne 0$, \eqref{d17z} yields
\begin{equation}
a_{k-3}=\frac{a_k-b_1a_{k-2}}{c_1}.\label{t5z}
\end{equation}
Using \eqref{icz}, \eqref{d10az} and \eqref{t5z}, we obtain
\begin{gather}
 a_{-1}=a_{-2}=0,\quad a_0=1,\label{b1z}\\
y_{-2}=y_{-1}=y_0=0,\quad y_1=1,\label{ic3z}\\
y_k=\sum_{j=0}^{k-1}a_j,\quad k\in\mathbb{N}.\label{ynz}
\end{gather}
By using \eqref{d17z}, \eqref{b1z}, \eqref{ynz} and Lemma \ref{lem2}, it follows
that closed form formulas for $a_k$ and $y_k$ can be found and consequently
a closed form formula for the solution to equation \eqref{r5a}.

Employing \eqref{d19z} in the second equation in \eqref{r1a} we obtain
\begin{equation}
w_n=\alpha^{dy_{n-2}}\beta^{1+bdy_{n-3}}z_{-1}^{da_{n-1}}z_0^{da_{n-2}}
w_0^{bda_{n-3}}.\label{b5}
\end{equation}
Some calculation show that formulas \eqref{d19z} and \eqref{b5} present
 a solution to system \eqref{r1a}, completing the proof of the theorem.
\end{proof}

The characteristic polynomial associated with equation \eqref{d17z} is 
\begin{equation}
p_3(\lambda)=\lambda^3-a\lambda-bd\,.\label{t27}
\end{equation}
The zeros of the polynomial are given by \eqref{s2g} with
\begin{equation}
\Delta_0=3a=:-3p\quad\text{and}\quad \Delta_1=27bd=:-27q.\label{de1}
\end{equation}
Hence, the analysis following Theorem \ref{thm3} also holds here,
 with the difference parameter $c$ replaced by $a$.
\smallskip

\noindent\textbf{Zeros of $p_3$ are different and none of them is 1.}
 In this case it must be $\Delta\ne 0$, that is,
$27(bd)^2\ne 4a^3$.
 If also $a+bd\ne 1$, then $p_3(1)\ne 0$. These conditions are satisfied, 
for example, if $a=bd\in\mathbb{N}$.
\smallskip

\noindent\textbf{Zeros of $p_3$ are different and one of them is 1.} 
Polynomial $p_3$ has a zero equal to 1 if $a+bd=1$. From this we have
\[
p_3(\lambda)=\lambda^3-a\lambda+a-1=(\lambda-1)(\lambda^2+\lambda+1-a),
\]
\begin{equation}
\lambda_{1,2}=\frac{-1\pm\sqrt{4a-3}}2,\quad \lambda_3=1,\label{l5z}
\end{equation}
are the zeros of \eqref{t27}. Since 1 is a simple zero, then it must
 be $p_3'(1)=3-a\ne 0$, that is, we also have $a\ne 3$.

\begin{corollary} \label{coro4} 
Assume that $a,b,d\in\mathbb{Z}$, $c=0$, $bd\ne0$, 
$\alpha, \beta, z_{-1}, z_0, w_0\in\mathbb{C}\setminus\{0\}$ and $\Delta\ne 0$. 
Then the following statements are true.
\begin{itemize}
\item[(a)]  If $a+bd\ne 1$, then the general solution to
\eqref{ms1} is given by \eqref{d19z} and \eqref{b5}, 
where $(a_n)_{n\ge-2}$ is given by
\eqref{s1m}, $(y_n)_{n\ge-2}$ is given by \eqref{a28b}, while 
$\lambda_j$-s, $j=\overline{1,3}$, are given by \eqref{s2g} and \eqref{de1}.

\item[(b)] If $a+bd=1$ and $a\ne 3$, then $p_3$ has a unique zero equal to 1, 
say $\lambda_3$, and the general solution to
\eqref{ms1} is given by formulas \eqref{d19z} and \eqref{b5},
 where $(a_n)_{n\ge-2}$ is given by
\eqref{s1m} with $\lambda_3=1$, $(y_n)_{n\ge-2}$ is given by \eqref{a29b},
 while $\lambda_j$-s, $j=\overline{1,3}$, are given by \eqref{l5z}.

\end{itemize}
\end{corollary}

\begin{corollary} \label{coro5} 
 Assume that $a, b, d\in\mathbb{Z}$, $c=0$, $bd\ne 0$,  
$\alpha, \beta, z_{-1}, z_0, w_0\in\mathbb{C}\setminus\{0\}$ and $\Delta=0$. 
Then the following statements are true.
\begin{itemize}
\item[(a)] If $a+bd\ne 1$, then the general solution to
\eqref{ms1} is given by \eqref{d19z} and \eqref{b5}, where 
$(a_n)_{n\ge-2}$ is given by
\eqref{a16}, $(y_n)_{n\ge-2}$ is given by \eqref{a32f}, while 
$\lambda_j$-s, $j=\overline{1,3}$, are given by \eqref{w7} and \eqref{de1}.

\item[(b)] If $a=3$ and $bd=-2$, then two zeros of \eqref{t27} are equal 
to 1, say, $\lambda_2$ and $\lambda_3$, and the general solution to system
\eqref{ms1} is given by \eqref{d19z} and \eqref{b5},  
where $(a_n)_{n\ge-2}$ is given by
\eqref{a16} with $\lambda_2=1$, $(y_n)_{n\ge-2}$ is given by \eqref{t29}, 
while $\lambda_1=-2$.

\item[(c)] Polynomial \eqref{t27} cannot have 1 as a simple zero.

\item[(d)] Polynomial \eqref{t27} cannot have a triple zero.

\end{itemize}
\end{corollary}

\begin{theorem} \label{thm5} 
 Assume that $a,b,c,d\in\mathbb{Z}$, $abcd\ne 0$, $\alpha, \beta, z_{-1}, 
z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. Then 
system \eqref{ms1} is solvable in closed form.
\end{theorem}

\begin{proof} 
The condition $\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$ 
along with \eqref{ms1} imply $z_nw_n\ne 0$ for $n\ge -1$. Hence
\begin{gather}
w_n^b=\frac{z_{n+1}}{\alpha z_{n-1}^a},\quad n\in\mathbb{N}_0,\label{d1} \\
w_{n+1}^b=\beta^b w_{n-1}^{bc}z_{n-1}^{bd},\quad
n\in\mathbb{N}_0.\label{d2}
\end{gather}
Using \eqref{d1} in \eqref{d2} we obtain
\begin{equation}
z_{n+2}=\alpha^{1-c}\beta^b z_n^{a+c}z_{n-1}^{bd}z_{n-2}^{-ac},\quad
n\in\mathbb{N}.\label{d3}
\end{equation}
Let $\delta=\alpha^{1-c}\beta^b$,
\begin{equation} a_1=0,\quad b_1=a+c,\quad c_1=bd,\quad d_1=-ac,\quad y_1=1.\label{ic}
\end{equation}
Then
\begin{equation}
z_{n+2}=\delta^{y_1}z_{n+1}^{a_1}z_n^{b_1}z_{n-1}^{c_1}z_{n-2}^{d_1},\quad
n\in\mathbb{N}.\label{d6}
\end{equation}
We have
\begin{equation}
\begin{aligned}
z_{n+2}
&=\delta^{y_1}(\delta z_n^{a_1}z_{n-1}^{b_1}z_{n-2}^{c_1}z_{n-3}^{d_1})^{a_1}
 z_n^{b_1}z_{n-1}^{c_1}z_{n-2}^{d_1},\\
&=\delta^{y_1+a_1}z_n^{a_1a_1+b_1}z_{n-1}^{b_1a_1+c_1}z_{n-2}^{c_1a_1+d_1}
 z_{n-3}^{d_1a_1}\\
&=\delta^{y_2}z_n^{a_2}z_{n-1}^{b_2}z_{n-2}^{c_2}z_{n-3}^{d_2},
\end{aligned} \label{d8}
\end{equation}
for $n\ge 2$, where
\begin{equation}
a_2:=a_1a_1+b_1,\; b_2:=b_1a_1+c_1,\; c_2:=c_1a_1+d_1,\; d_2:=d_1a_1,\;
y_2:=y_1+a_1.\label{d85}
\end{equation}
Assume
\begin{equation}
z_{n+2}=\delta^{y_k}z_{n+2-k}^{a_k}z_{n+1-k}^{b_k}z_{n-k}^{c_k}z_{n-k-1}^{d_k},\label{d9}
\end{equation}
for a $k\ge 2$ and every $n\ge k$, and that
\begin{gather} \label{d10}
\begin{gathered}
a_k=a_1a_{k-1}+b_{k-1},\quad b_k=b_1a_{k-1}+c_{k-1},\\
c_k=c_1a_{k-1}+d_{k-1},\quad d_k=d_1a_{k-1},
\end{gathered} \\
 y_k=y_{k-1}+a_{k-1}.\label{d10a}
\end{gather}
Using \eqref{d6} in \eqref{d9}, it follows that
\begin{equation}
\begin{aligned}
z_{n+2}
&=\delta^{y_k}(\delta z_{n+1-k}^{a_1}z_{n-k}^{b_1}z_{n-k-1}^{c_1}
 z_{n-k-2}^{d_1})^{a_k}z_{n+1-k}^{b_k}z_{n-k}^{c_k}z_{n-k-1}^{d_k}\\
&=\delta^{y_k+a_k}z_{n+1-k}^{a_1a_k+b_k}z_{n-k}^{b_1a_k+c_k}z_{n-k-1}^{c_1a_k+d_k}
 z_{n-k-2}^{d_1a_k}\\
&=\delta^{y_{k+1}}z_{n+1-k}^{a_{k+1}}z_{n-k}^{b_{k+1}}z_{n-k-1}^{c_{k+1}}
 z_{n-k-2}^{d_{k+1}},
\end{aligned} \label{d11}
\end{equation}
for $n\ge k+1$, where
\begin{gather} \label{d12}
\begin{gathered}
a_{k+1}:=a_1a_k+b_k,\quad b_{k+1}:=b_1a_k+c_k,\\
c_{k+1}:=c_1a_k+d_k,\quad d_{k+1}:=d_1a_k,
\end{gathered} \\
 y_{k+1}:=y_k+a_k.\label{d12t}
\end{gather}
From \eqref{d8}, \eqref{d85}, \eqref{d11}-\eqref{d12t}, the inductive
argument shows that \eqref{d9}-\eqref{d10a} hold for
$2\le k\le n$.

Setting $k=n$ in \eqref{d9}, then using \eqref{d10} and \eqref{d10a} 
in the obtained equality, we have
\begin{equation}
\begin{aligned}
z_{n+2}
=&\delta^{y_n}z_2^{a_n}z_1^{b_n}z_0^{c_n}z_{-1}^{d_n}\\
=&(\alpha^{1-c}\beta^b)^{y_n}(\alpha\beta^bw_{-1}^{bc}z_{-1}^{bd}z_0^a)^{a_n}(\alpha z_{-1}^a
 w_0^b)^{b_n}z_0^{c_n}z_{-1}^{d_n}\\
=&\alpha^{(1-c)y_n+a_n+b_n}\beta^{by_n+ba_n} z_{-1}^{bda_n+ab_n+d_n}z_0^{aa_n+c_n}
 w_{-1}^{bca_n}w_0^{bb_n}\\
=&\alpha^{y_{n+2}-cy_n}\beta^{by_{n+1}}z_{-1}^{a_{n+3}-ca_{n+1}}z_0^{a_{n+2}-ca_n}
 w_{-1}^{bca_n}w_0^{ba_{n+1}},\quad
n\in\mathbb{N}.
\end{aligned} \label{d19}
\end{equation}
From \eqref{d10} we obtain
 \begin{equation}
a_k=b_1a_{k-2}+c_1a_{k-3}+d_1a_{k-4},\quad k\geq 5 .\label{d17}
\end{equation}
Equalities \eqref{d10a} and \eqref{d17} yield
\begin{gather}\label{d50}
a_{-3}=a_{-2}=a_{-1}=0,\quad a_0=1; \\
y_{-3}=y_{-2}=y_{-1}=y_0=0,\quad y_1=1,\label{ic3} \\
y_k=\sum_{j=0}^{k-1}a_j.\label{yk}
\end{gather}

Since the initial-value problem \eqref{d17}-\eqref{d50} is solvable, 
a formula for $a_k$ can be found. Using it in \eqref{yk} and applying 
Lemma \ref{lem2},  a formula for $y_k$ is found. By using these 
two formulas in \eqref{d19} 
we obtain a formula for $z_n$. So, \eqref{d3} is solvable.


On the other hand, we have
\begin{gather} 
z_{n-1}^d=\frac{w_{n+1}}{\beta w_{n-1}^c},\quad
n\in\mathbb{N}_0,\label{f1}\\
 z_{n+1}^d=\alpha^dz_{n-1}^{ad}w_n^{bd},\quad
n\in\mathbb{N}_0,\label{f2}
\end{gather}
and consequently
\begin{equation}
w_{n+3}=\alpha^d\beta^{1-a}w_{n+1}^{a+c}w_n^{bd}w_{n-1}^{-ac},\quad
n\in\mathbb{N}_0.\label{f3}
\end{equation}
As above, we have
\begin{equation}
w_{n+3}=\eta^{y_k}w_{n+3-k}^{a_k}w_{n+2-k}^{b_k}w_{n+1-k}^{c_k}w_{n-k}^{d_k},\quad
n\ge k-1,\label{f9}
\end{equation}
where $\eta=\alpha^d\beta^{1-a}$, $(a_k)_{k\in\mathbb{N}}$,
$(b_k)_{k\in\mathbb{N}}$, $(c_k)_{k\in\mathbb{N}}$ and $(d_k)_{k\in\mathbb{N}}$
 are defined by \eqref{ic} and
\eqref{d10}, while $(y_k)_{k\in\mathbb{N}}$ is defined by \eqref{d10a}
and \eqref{ic3}.

From \eqref{f9} with $k=n+1$ and by using \eqref{z}, we obtain
\begin{equation}
\begin{aligned}
w_{n+3}
&=\eta^{y_{n+1}}w_2^{a_{n+1}}w_1^{b_{n+1}}w_0^{c_{n+1}}w_{-1}^{d_{n+1}}\\
&=(\alpha^d\beta^{1-a})^{y_{n+1}}(\beta w_0^cz_0^d)^{a_{n+1}}(\beta w_{-1}^cz_{-1}^d)^{b_{n+1}}
 w_0^{c_{n+1}}w_{-1}^{d_{n+1}}\\
&=\alpha^{dy_{n+1}}\beta^{(1-a)y_{n+1}+a_{n+1}+b_{n+1}}z_{-1}^{db_{n+1}}z_0^{da_{n+1}}
 w_{-1}^{cb_{n+1}+d_{n+1}}w_0^{ca_{n+1}+c_{n+1}}\\
&=\alpha^{dy_{n+1}}\beta^{y_{n+3}-ay_{n+1}}z_{-1}^{da_{n+2}}z_0^{da_{n+1}}w_{-1}
 ^{c(a_{n+2}-aa_n)}w_0^{a_{n+3}-aa_{n+1}},\quad
\end{aligned} \label{f19}
\end{equation}
for $n\in\mathbb{N}_0.$

The formulas for $a_k$ and $y_k$ are obtained as above. Using them in 
\eqref{f19} is get a formula for a solution to \eqref{f3}. By some calculation 
it is checked that \eqref{d19} and \eqref{f19} are formulas for a solution 
to \eqref{ms1}. Thus,
the system is also solvable, as claimed.
\end{proof}



\begin{corollary} \label{coro6}
 Assume that $a,b,c,d\in\mathbb{Z}$, $abcd\ne 0$, 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$.  
Then  the general solution to system
\eqref{ms1} is given by \eqref{d19} and \eqref{f19}, where 
$(a_k)_{k\in\mathbb{N}}$ is defined by \eqref{d17} and \eqref{d50}, 
while $(y_k)_{k\in\mathbb{N}}$ is defined by \eqref{ic3} and \eqref{yk}.
\end{corollary}

The characteristic polynomial associated to \eqref{d17} is
\begin{equation}
p_4(\lambda)=\lambda^4-(a+c)\lambda^2-bd\lambda+ac.\label{q28}
\end{equation}
Since $ac\ne 0$, it is of the forth order.

The equation $p_4(\lambda)=0$ can be written as
\begin{equation}
\Big(\lambda^2-\frac{a+c}2+\frac{s}2\Big)^2-\Big(s\lambda^2+bd\lambda
+\Big(\frac{s-(a+c)}2\Big)^2-ac\Big)=0.\label{fe}
\end{equation}
Let $s$ satisfy $(bd)^2=s(s-a-c)^2-4acs$, that is,
\begin{equation}
s^3-2(a+c)s^2+(a-c)^2s-(bd)^2=0.\label{te}
\end{equation}
For such $s$, \eqref{fe} becomes
$$
\Big(\lambda^2-\frac{a+c}2+\frac{s}2\Big)^2
-\Big(\sqrt{s}\;\lambda+\frac{bd}{2\sqrt{s}}\Big)^2=0,
$$
which is equivalent to the following two quadratic equations
\begin{gather}
\lambda^2-\sqrt{s}\,\lambda+\frac{s-a-c}2-\frac{bd}{2\sqrt{s}}=0,\label{feq1}\\
\lambda^2+\sqrt{s}\,\lambda+\frac{s-a-c}2+\frac{bd}{2\sqrt{s}}=0.\label{feq2}
\end{gather}
By using the change of variables $s=t+\frac{2(a+c)}3$ in \eqref{te} we obtain
\begin{equation}
t^3-\Big(\frac{4(a+c)^2}3-(a-c)^2\Big)t-\frac{16(a+c)^3}{27}
+\frac{2(a+c)(a-c)^2}3-(bd)^2=0.\label{to}
\end{equation}
Let
$$
p=-\frac{a^2+14ac+c^2}3\quad\text{and}\quad
q=\frac{2a^3+2c^3-66ac(a+c)-27(bd)^2}{27}.
$$
Then, we can choose $t$ as one of the three possible values of the quantity
\begin{equation}
t=\sqrt[3]{-\frac{q}2-\sqrt{\frac{q^2}4+\frac{p^3}{27}}}
 +\sqrt[3]{-\frac{q}2+\sqrt{\frac{q^2}4+\frac{p^3}{27}}}.\label{te2}
\end{equation}
If we use the change $p=-\Delta_0/3$ and $q=-\Delta_1/27$ in \eqref{te2} we obtain
$$
t=\frac{1}{3\sqrt[3]{2}}\Big(\sqrt[3]{\Delta_1-\sqrt{\Delta_1^2-4\Delta_0^3}}
+\sqrt[3]{\Delta_1+\sqrt{\Delta_1^2-4\Delta_0^3}}\Big).
$$
Solutions to \eqref{feq1} and \eqref{feq2} are the zeros of polynomial \eqref{q28}
 and they are
\begin{gather}
\lambda_1=\frac{1}{2}\sqrt{t+\frac{2(a+c)}3}
+\frac{1}{2}\sqrt{\frac{4(a+c)}3-t-\frac{Q}{4\sqrt{t+\frac{2(a+c)}3}}},\label{l11}
\\
\lambda_2=\frac{1}{2}\sqrt{t+\frac{2(a+c)}3}
-\frac{1}{2}\sqrt{\frac{4(a+c)}3-t-\frac{Q}{4\sqrt{t+\frac{2(a+c)}3}}},\label{l12}
\\
\lambda_3=-\frac{1}{2}\sqrt{t+\frac{2(a+c)}3}
+\frac{1}{2}\sqrt{\frac{4(a+c)}3-t+\frac{Q}{4\sqrt{t+\frac{2(a+c)}3}}},\label{l13}
\\
\lambda_4=-\frac{1}{2}\sqrt{t+\frac{2(a+c)}3}
-\frac{1}{2}\sqrt{\frac{4(a+c)}3-t+\frac{Q}{4\sqrt{t+\frac{2(a+c)}3}}},\label{l14}
\end{gather}
where
\begin{equation}
Q:=-8bd.\label{l21}
\end{equation}

The nature of $\lambda_j$, $j=\overline{1,4}$, depends on the sign of
\begin{equation} \Delta:=\frac{1}{27}(4\Delta_0^3-\Delta_1^2),\label{l22}\end{equation}
where
\begin{align}
\Delta_0:=&a^2+14ac+c^2,\label{l21a}\\
\Delta_1:=&-2a^3-2c^3+66ac(a+c)+27(bd)^2,\label{l21b}
\end{align}
and the signs of
\begin{equation} P:=-8(a+c)\label{23a}\end{equation}
and
\begin{equation} D:=-16(a-c)^2.\label{23}
\end{equation}
\smallskip

\noindent\textbf{Zeros of $p_4$ are different and none of them is 1.} 
If $a$, $b$, $c$ and $d$ are chosen such that $\Delta_0<0$, that is, 
$a^2+14ac+c^2<0$, then it will be $\Delta<0$, from which by Lemma 
\ref{lem3} it follows that $p_4$ has four different zeros. Let $t_{1,2}=-7\pm\sqrt{48}$. 
Then if we choose $a,c\in\mathbb{Z}\setminus\{0\}$ such that 
$a/c\in(t_1,t_2)$, then we have such a situation.
\smallskip

\noindent\textbf{Zeros of $p_4$ are different and one of them is 1.} 
 Polynomial \eqref{q28} has a zero equal to 1 if
$p_4(1)=1-a-c-bd+ac=0$, that is, if 
\begin{equation}
(a-1)(c-1)=bd,\label{w1}
\end{equation}
so that
\begin{equation}
p_4(\lambda)=\lambda^4-(a+c)\lambda^2-(a-1)(c-1)\lambda +ac.
\end{equation}
Thus, if we choose $a$ and $c$ such that $p_4'(1)=3-a-c-ac\ne 0$,
that is, $(a+1)(c+1)\ne 4$, then $p_4$ will be such a polynomial
if $\Delta\ne 0$. For example, if $a=-3$ and $c=2$, then
$bd=-4\ne 0$, $\Delta\ne 0$ which means that the characteristic polynomial has
all zeros mutually different and exactly one of them is equal to 1
\begin{equation}
p_4(\lambda)=\lambda^4+\lambda^2+4\lambda-6=(\lambda-1)(\lambda^3+\lambda^2+2\lambda+6).
\end{equation}

Since in these two cases $\lambda_j\ne \lambda_i$, $i\ne j$, then the general 
solution to \eqref{d17} is
\begin{equation}
a_n=\gamma_1\lambda_1^n+\gamma_2\lambda_2^n+\gamma_3\lambda_3^n+\gamma_4\lambda_4^n,\quad
n\in\mathbb{N},\label{d51}
\end{equation}
where $\gamma_i$, $i=\overline{1,4}$, are arbitrary constants.

The solution to \eqref{d17} satisfying \eqref{d50} is
\begin{equation}
\begin{aligned}
a_n=&\sum_{j=1}^4\frac{\lambda_j^{n+3}}{p_4'(\lambda_j)}\\
&=\frac{\lambda_1^{n+3}}{(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)(\lambda_1-\lambda_4)}
 +\frac{\lambda_2^{n+3}}{(\lambda_2-\lambda_1)(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)}\\
&\quad+\frac{\lambda_3^{n+3}}{(\lambda_3-\lambda_1)(\lambda_3-\lambda_2)(\lambda_3-\lambda_4)}
+\frac{\lambda_4^{n+3}}{(\lambda_4-\lambda_1)(\lambda_4-\lambda_2)(\lambda_4-\lambda_3)},
\end{aligned} \label{s1}
\end{equation}
for $n\ge -3$ (\cite{jia2015-pts}).

From \eqref{yk} and \eqref{s1} it follows that
\begin{equation}
y_n=\sum_{j=0}^{n-1}\sum_{i=1}^4\frac{\lambda_i^{j+3}}{p_4'(\lambda_i)}
=\sum_{i=1}^4\frac{\lambda_i^3(\lambda_i^n-1)}{p_4'(\lambda_i)(\lambda_i-1)},
\quad n\in\mathbb{N},\label{yn}
\end{equation}
when $\lambda_i\ne 1$, $i=\overline{1,4}$, and
\begin{equation}
y_n=\frac{n}{3-a-c-ac}+\sum_{i=2}^4
\frac{\lambda_i^3(\lambda_i^n-1)}{p_4'(\lambda_i)(\lambda_i-1)},\quad n\in\mathbb{N},\label{yn1}
\end{equation}
when one of the zeros of $p_4$ is 1 (here $\lambda_1=1$).

Note that if one of the zeros of $p_4$ is 1, then we have
\begin{equation}
p_4(\lambda)=(\lambda-1)(\lambda^3+\lambda^2-(a+c-1)\lambda-ac).\label{w6}
\end{equation}
The change of variables $\lambda=t-\frac13$ transforms the following equation
$$
\lambda^3+\lambda^2-(a+c-1)\lambda-ac=0,
$$
into
\begin{equation}
t^3+\tilde pt+\tilde q=0,\label{w8}
\end{equation}
where
\begin{equation}
\tilde p=\frac23-a-c\quad\text{and}\quad \tilde q=\frac{9(a+c)-7-27ac}{27}.\label{pq}
\end{equation}
The zeros of \eqref{w8} are
$$
t_l=\varepsilon^{l-1}\sqrt[3]{-\frac{\tilde q}2
-\sqrt{\frac{\tilde q^2}4+\frac{\tilde p^3}{27}}}
+\varepsilon^{l-1}\sqrt[3]{-\frac{\tilde q}2
+\sqrt{\frac{\tilde q^2}4+\frac{\tilde p^3}{27}}},\quad l=\overline{1,3},
$$
where $\varepsilon^3=1$ and $\varepsilon \ne 1$.
Hence,
\begin{equation}
\lambda_j=-\frac{1}3+\varepsilon^{j-2}\sqrt[3]{-\frac{\tilde q}2
-\sqrt{\frac{\tilde q^2}4+\frac{\tilde p^3}{27}}}
+\varepsilon^{j-2}\sqrt[3]{-\frac{\tilde q}2
+\sqrt{\frac{\tilde q^2}4+\frac{\tilde p^3}{27}}},\label{l1}
\end{equation}
for $j=\overline{2,4}$, are the other three zeros of $p_4$, in this case.

The previous analysis along with Corollary \ref{coro6} implies 
the following corollary.

\begin{corollary} \label{coro7}
 Assume that $a,b,c,d\in\mathbb{Z}$, $abcd\ne 0$, 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$ and 
$\Delta\ne 0$. Then the following statements are true.
\begin{itemize}
\item[(a)]  If $(a-1)(c-1)\ne bd$, then the general solution to 
\eqref{ms1} is given by \eqref{d19} and \eqref{f19}, where 
$(a_n)_{n\ge-3}$ is given by \eqref{s1}, $(y_n)_{n\ge -3}$ 
is given by \eqref{yn}, while $\lambda_j$-s, $j=\overline{1,4}$, are given by
 \eqref{l11}-\eqref{l14}.

\item[(b)] If $(a-1)(c-1)=bd$ and $(a+1)(c+1)\ne 4$, then the general 
solution to \eqref{ms1} is given by \eqref{d19} and \eqref{f19}, 
where $(a_n)_{n\ge-3}$ is given by \eqref{s1} with $\lambda_1=1$,
 $(y_n)_{n\ge -3}$ is given by \eqref{yn1}, $\lambda_1=1$, while 
$\lambda_j$-s, $j=\overline{2,4}$, are given by \eqref{l1} and \eqref{pq}.

\end{itemize}
\end{corollary} 
\smallskip

\noindent\textbf{$p_4$ has only one double zero which is equal to 1.} 
Polynomial \eqref{q28} will have a double zero equal to 1 if \eqref{w1} holds and
\begin{equation}
(a+1)(c+1)=4.\label{w3}
\end{equation}
From \eqref{w3} we have that one of the following cases must occur:
(1) $a=3$ and $c=0$, (2) $a=0$ and $c=3$; (3) $a=c=1$; (4) $a=-5$ and $c=-2$;
(5) $a=-2$ and $c=-5$; (6) $a=c=-3$.
If $a=0$ or $c=0$, then $ac=0$, which is a contradiction.

If $a=c=1$, then
$$p_4(\lambda)=\lambda^4-2\lambda^2+1=(\lambda-1)^2(\lambda+1)^2,$$
and
\begin{equation}
\lambda_{1,2}=1,\quad \lambda_{3,4}=-1.\label{c11}
\end{equation}

If $a=-5$ and $c=-2$ or $a=-2$ and $c=-5$, then
$$
p_4(\lambda)=\lambda^4+7\lambda^2-18\lambda^2+10=(\lambda-1)^2(\lambda^2+2\lambda+10),
$$
and
\begin{equation}
\lambda_{1,2}=1,\quad \lambda_{3,4}=-1\pm 3i.\label{c12}
\end{equation}

If $a=c=-3$, then
$$
p_4(\lambda)=\lambda^4+6\lambda^2-16\lambda+9=(\lambda-1)^2(\lambda^2+2\lambda+9),
$$
and
\begin{equation}
\lambda_{1,2}=1,\quad \lambda_{3,4}=-1\pm2\sqrt2 i.\label{c14}
\end{equation}

From this, we have proved in passing, that there are no such
 $a,c\in\mathbb{Z}\setminus\{0\}$, such that 1 is a triple zero of $p_4$, 
or that $p_4$ has two pairs of double zeros such that one of them is equal to 1.

In these four cases we have (see, for example, \cite{ejqtde-4th-2017})
\begin{equation}
\begin{aligned}
a_n&=\frac{n(1-\lambda_3)(1-\lambda_4)+3\lambda_3\lambda_4-2\lambda_3-2\lambda_4+1}{(1-\lambda_3)^2(1-\lambda_4)^2}\\
&\quad +\frac{\lambda_3^{n+3}}{(\lambda_3-1)^2(\lambda_3-\lambda_4)}
 +\frac{\lambda_4^{n+3}}{(\lambda_4-1)^2(\lambda_4-\lambda_3)},
\end{aligned} \label{an3}
\end{equation}
and
\begin{equation}
\begin{aligned}
y_n&=\sum_{j=0}^{n-1}\Big(\frac{j(1-\lambda_3)(1-\lambda_4)+3\lambda_3\lambda_4-2\lambda_3
 -2\lambda_4+1}{(1-\lambda_3)^2(1-\lambda_4)^2}\\
&\quad+\frac{\lambda_3^{j+3}}{(\lambda_3-1)^2(\lambda_3-\lambda_4)}
 +\frac{\lambda_4^{j+3}}{(\lambda_4-1)^2(\lambda_4-\lambda_3)}\Big)\\
&=\frac{(n-1)n}{2(1-\lambda_3)(1-\lambda_4)}
 +\frac{n(3\lambda_3\lambda_4-2\lambda_3-2\lambda_4+1)}{(1-\lambda_3)^2(1-\lambda_4)^2}\\
&\quad +\frac{\lambda_3^3(\lambda_3^n-1)}{(\lambda_3-1)^3(\lambda_3-\lambda_4)}
 +\frac{\lambda_4^3(\lambda_4^n-1)}{(\lambda_4-1)^3(\lambda_4-\lambda_3)}.
\end{aligned} \label{yn3}
\end{equation}
\smallskip


\noindent\textbf{Exactly one double zero which is different from 1.}
 According to Lemma \ref{lem3}, in this case it must be $\Delta=0$, that is,
\begin{equation}
(a^2+14ac+c^2)^3=\Big(a^3+c^3-33ac(a+c)-\frac{27}2(bd)^2\Big)^2,\label{c15}
\end{equation}
and that
\begin{equation}
(a-1)(c-1)\ne bd,\quad a+c>0,\quad a\ne c.\label{c16}
\end{equation}

The problem of the existence $a, c\in\mathbb{Z}\setminus\{0\}$ such that 
\eqref{c15} and \eqref{c16} hold seems quite technical and we leave it to 
the reader as an open problem.

Since, in the case $\lambda_1=\lambda_2$, $\lambda_i\ne \lambda_j$, $2\le i,j\le 4$, we have 
that the general solution to \eqref{d17} has the form
 \begin{equation}
a_n=(\gamma_1+\gamma_2n)\lambda_2^n+\gamma_3\lambda_3^n+\gamma_4\lambda_4^n, \quad
n\in\mathbb{N},\label{c2}
\end{equation}
where $\gamma_i$, $i=\overline{1,4}$, are arbitrary constants, and the solution
satisfying \eqref{d50} can be obtained, for example, by letting $\lambda_1\to\lambda_2$
in \eqref{s1} \cite{ejqtde-4th-2017},
\begin{equation}
\begin{aligned}
a_n&=\frac{\lambda_2^{n+2}((n+3)(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)-\lambda_2(2\lambda_2-\lambda_3-\lambda_4))
}{(\lambda_2-\lambda_3)^2(\lambda_2-\lambda_4)^2}\\
&\quad +\frac{\lambda_3^{n+3}}{(\lambda_3-\lambda_2)^2(\lambda_3-\lambda_4)}
 +\frac{\lambda_4^{n+3}}{(\lambda_4-\lambda_2)^2(\lambda_4-\lambda_3)}.
\end{aligned} \label{s1a}
\end{equation}
From \eqref{yk}, \eqref{s1a} and by Lemma \ref{lem2}, we obtain
\begin{equation}
\begin{aligned}
y_n&=\sum_{j=0}^{n-1}\Big(\frac{\lambda_2^{j+2}((j+3)(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)
 -\lambda_2(2\lambda_2-\lambda_3-\lambda_4))}{(\lambda_2-\lambda_3)^2(\lambda_2-\lambda_4)^2}\\
&\quad+\frac{\lambda_3^{j+3}}{(\lambda_3-\lambda_2)^2(\lambda_3-\lambda_4)}
 +\frac{\lambda_4^{j+3}}{(\lambda_4-\lambda_2)^2(\lambda_4-\lambda_3)}\Big)\\
&=\frac{\lambda_2^3-n\lambda_2^{n+2}+(n-1)\lambda_2^{n+3}}{(\lambda_2-\lambda_3)(\lambda_2-\lambda_4)(1-\lambda_2)^2}
+\frac{(\lambda_2^4-2\lambda_2^3\lambda_3-2\lambda_2^3\lambda_4+3\lambda_2^2\lambda_3\lambda_4)(\lambda_2^n-1)}
 {(\lambda_2-\lambda_3)^2(\lambda_2-\lambda_4)^2(\lambda_2-1)}\\
&\quad+\frac{\lambda_3^3(\lambda_3^n-1)}{(\lambda_3-\lambda_2)^2(\lambda_3-\lambda_4)(\lambda_3-1)}
 +\frac{\lambda_4^3(\lambda_4^n-1)}{(\lambda_4-\lambda_2)^2(\lambda_4-\lambda_3)(\lambda_4-1)}.
\end{aligned}\label{yn2}
\end{equation}

From the previous analysis and Corollary \ref{coro6} we obtain the following result.

\begin{corollary} \label{coro8}
 Assume that $a,b,c,d\in\mathbb{Z}$, $abcd\ne 0$ and 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. 
Then the following statements are true.
\begin{itemize}
\item[(a)] If only one of the zeros of polynomial \eqref{q28} is double 
and different from $1$, then the general solution to
\eqref{ms1} is given by \eqref{d19} and \eqref{f19}, where 
$(a_n)_{n\ge-3}$ is given by \eqref{s1a}, while $(y_n)_{n\ge -3}$ 
is given by \eqref{yn2}.

\item[(b)] If $1$ is a unique double zero of polynomial $p_4$, say 
$\lambda_1=\lambda_2=1$, then the general solution to
\eqref{ms1} is given by \eqref{d19} and \eqref{f19}, where 
$(a_n)_{n\ge-3}$ is given by \eqref{an3}, $(y_n)_{n\ge -3}$ is given 
by \eqref{yn3}, while $\lambda_{3,4}$ are given by \eqref{c11} if $a=c=1$, 
by \eqref{c12} if $a=-5$, $c=-2$ or $a=-2$, $c=-5$, and by 
\eqref{c14} if $a=c=-3$.

\end{itemize}
\end{corollary}
\smallskip

\noindent\textbf{Two pairs of different double zeros.}
 Let $bd=0$, then
\begin{gather*}
p_4(\lambda)=\lambda^4-(a+c)\lambda^2+ac=(\lambda^2-a)(\lambda^2-c),\\
\lambda_{1,2}=\pm\sqrt{a},\quad \lambda_{3,4}=\pm\sqrt{c}.
\end{gather*}
Hence, if $a=c\ne 0$, $p_4$ in this case has two pairs of different double zeros
\begin{equation}
\lambda_{1,3}=\sqrt{a},\quad \lambda_{3,4}=-\sqrt{a}.\label{c18}
\end{equation}
The general solution to \eqref{d17} in this case has the  form
\begin{equation}
a_n=(\gamma_1+\gamma_2n)\lambda_2^n+(\gamma_3+\gamma_4n)\lambda_4^n, \quad n\in\mathbb{N},\label{d53}
\end{equation}
where $\gamma_i$, $i=\overline{1,4}$ are constants.

The solution to \eqref{d17} of the form in \eqref{d53} and satisfying \eqref{d50}
 is given by \cite{ejqtde-4th-2017}
\begin{equation}
\begin{aligned}
a_n&=\frac{\lambda_2^{n+2}(n(\lambda_2-\lambda_4)^2+\lambda_2^2-4\lambda_2\lambda_4+3\lambda_4^2)
}{(\lambda_2-\lambda_4)^4}\\
&\quad +\frac{\lambda_4^{n+2}(n(\lambda_4-\lambda_2)^2+\lambda_4^2-4\lambda_2\lambda_4
+3\lambda_2^2)}{(\lambda_4-\lambda_2)^4}.
\end{aligned}\label{an1}
\end{equation}
From \eqref{yk}, \eqref{an1} and by Lemma \ref{lem2}, we obtain
\begin{equation}
\begin{aligned}
y_n
&=\sum_{j=0}^{n-1}\Big(\frac{\lambda_2^{j+2}(j(\lambda_2-\lambda_4)^2+\lambda_2^2-4\lambda_2\lambda_4
 +3\lambda_4^2)}{(\lambda_2-\lambda_4)^4}\\
&\quad+\frac{\lambda_4^{j+2}(j(\lambda_4-\lambda_2)^2+\lambda_4^2-4\lambda_2\lambda_4+3\lambda_2^2)}
 {(\lambda_4-\lambda_2)^4}\Big)\\
&=\frac{\lambda_2^3-n\lambda_2^{n+2}+(n-1)\lambda_2^{n+3}}{(\lambda_2-\lambda_4)^2(1-\lambda_2)^2}
+\frac{(\lambda_2^4-4\lambda_2^3\lambda_4+3\lambda_2^2\lambda_4^2)(\lambda_2^n-1)}{(\lambda_2-\lambda_4)^4(\lambda_2-1)}\\
&\quad +\frac{\lambda_4^3-n\lambda_4^{n+2}+(n-1)\lambda_4^{n+3}}{(\lambda_4-\lambda_2)^2(1-\lambda_4)^2}
 + \frac{(\lambda_4^4-4\lambda_2\lambda_4^3+3\lambda_2^2\lambda_4^2)(\lambda_4^n-1)}{(\lambda_4-\lambda_2)^4(\lambda_4-1)}.
\end{aligned} \label{yn5}
\end{equation}

\begin{corollary} \label{coro9} 
 Assume that $a,b,c,d\in\mathbb{Z}$, $ac\ne 0$ and 
$\alpha, \beta, z_{-1}, z_0, w_{-1}, w_0\in\mathbb{C}\setminus\{0\}$. 
Then the following statements are true.
\begin{itemize}
\item[(a)] If polynomial $p_4$ has two pairs of double zeros both different 
from $1$, then the general solution to
\eqref{ms1} is given by \eqref{d19} and \eqref{f19}, where
 $(a_n)_{n\ge-3}$ is given by \eqref{an1}, while $(y_n)_{n\ge -3}$
 is given by \eqref{yn5}. If additionally  $bd=0$ and $a=c$, then 
$\lambda_j$-s, $j=\overline{1,4}$, are given by \eqref{c18}.

\item[(b)] For $a=c=1$ and $bd=0$ the polynomial in \eqref{q28} can have 
two pairs of double zeros such that one of them is equal to $1$, and the 
general solution to \eqref{ms1} is given by \eqref{d19} and \eqref{f19}, 
where $(a_n)_{n\ge-3}$ is given by \eqref{an1}, while $(y_n)_{n\ge -3}$ 
is given by \eqref{yn5} with $\lambda_2=1$ and $\lambda_4=-1$.

\end{itemize}
\end{corollary}
\smallskip


\noindent\textbf{Triple zero case.} 
Since $\Delta_0\ne 0$, when $a, c\in\mathbb{Z}\setminus\{0\}$, by 
Lemma \ref{lem3} it follows that $p_4$ cannot have a triple zero. 
Consequently, it cannot have a quadruple zero.


\begin{thebibliography}{00}

\bibitem{ra} R.~P.~Agarwal; 
\emph{Difference Equations and Inequalities: Theory, Methods, and Applications}
 2nd Edition, Marcel Dekker Inc., New York, Basel, 2000.

\bibitem{am1} A.~Andruch-Sobilo, M.~Migda;
Further properties of the rational recursive sequence 
$x_{n+1}=ax_{n-1}/(b+cx_nx_{n-1})$, \emph{Opuscula Math.} \textbf{26} (3) (2006),
387-394.

\bibitem{amc218-sde} L.~Berg, S.~Stevi\' c;
On some systems of difference equations, \emph{Appl. Math. Comput.} 
\textbf{218} (2011), 1713-1718.

\bibitem{f} D.~K.~Faddeyev;
\emph{Lectures on Algebra}, Nauka, Moscow, 1984 (in Russian).

\bibitem{fp1} N.~Fotiades, G.~Papaschinopoulos;
Asymptotic behavior of the positive solutions of a system of $k$ difference 
equations of exponential form; \emph{Dynam. Contin. Discrete Impuls. Systems
Ser. A} \textbf{19} (5) (2012), 585-597.

\bibitem{i5} B.~Iri\v{c}anin, S.~Stevi\' c;
 Eventually constant solutions of a rational difference equation, 
\emph{Appl. Math. Comput.}  \textbf{215} (2009), 854-856.

\bibitem{j} C.~Jordan; 
\emph{Calculus of Finite Differences}, Chelsea Publishing Company,
 New York, 1956.

\bibitem{kk} G.~L.~Karakostas, Convergence of a difference equation via the full limiting sequences method; \emph{Differ. Equ. Dyn. Syst.} \textbf{1} (4) (1993), pp. 289-294.

\bibitem{k} V.~A.~Krechmar; 
\emph{A Problem Book in Algebra}, Mir Publishers, Moscow, 1974.

\bibitem{ll} H.~Levy, F.~Lessman; 
\emph{Finite Difference Equations}, Dover Publications, New York, 1992.

\bibitem{jm} J.~Migda;
 Qualitative approximation of solutions to difference equations,
\emph{Electron. J. Qual. Theory Differ. Equ.} Vol. 2015, Article no. 32, (2015), 
26 pages.

\bibitem{mk} D.~S.~Mitrinovi\'c, J.~D.~Ke\v{c}ki\'c;
\emph{Methods for Calculating Finite Sums}, Nau\v{c}na Knjiga, Beograd,
1984 (in Serbian).

\bibitem{p3-ejqtde} G.~Papaschinopoulos, N.~Psarros, K.~B.~Papadopoulos;
 On a system of $m$ difference equations having exponential terms,
 \emph{Electron. J. Qual. Theory Differ. Equ.} Vol. 2015, Article no. 5,
(2015), 13 pages.

\bibitem{ps1} G.~Papaschinopoulos, C.~J.~Schinas;
 On a system of two nonlinear difference equations,
 \emph{J. Math. Anal. Appl.} \textbf{219} (2) (1998), 415-426.

\bibitem{ps2} G.~Papaschinopoulos, C.~J.~Schinas;
 On the behavior of the solutions of a system of two nonlinear difference equations; 
\emph{Comm. Appl. Nonlinear Anal.} \textbf{5} (2) (1998), 47-59.

\bibitem{ps3} G.~Papaschinopoulos, C.~J.~Schinas;
 Invariants for systems of two nonlinear difference
equations, \emph{Differential Equations Dynam. Systems} \textbf{7} (2)
(1999), 181-196.

\bibitem{ps4} G.~Papaschinopoulos, C.~J.~Schinas;
 Invariants and oscillation for systems of two nonlinear difference
equations, \emph{Nonlinear Anal. TMA} \textbf{46} (7) (2001), 967-978.

\bibitem{ps-jdea2011} G.~Papaschinopoulos, C.~J.~Schinas;
 Oscillation and asymptotic stability of two systems of
difference equations of rational form,
 \emph{J. Difference Equat. Appl.} \textbf{7} (2001), 601-617.

\bibitem{psh} G.~Papaschinopoulos, C.~Schinas, V.~Hatzifilippidis;
 Global behavior of the solutions of a max-equation and of a system
of two max-equations, \emph{J. Comput. Anal. Appl.} \textbf{5} (2)
(2003), 237-254.

\bibitem{pst2} G. Papaschinopoulos, G. Stefanidou;
 Asymptotic behavior of the solutions of a class of rational difference equations,
 \emph{Inter. J. Difference Equations} \textbf{5} (2) (2010), 233-249.

\bibitem{rv} V.~D.~R\u{a}dulescu;
 Nonlinear elliptic equations with variable exponent: old and new,
\emph{Nonlinear Anal.} \textbf{121} (2015), 336-369.


\bibitem{r} E.~L.~Rees;
 Graphical discussion of the roots of a quartic equation,
 \emph{Amer. Math. Monthly} \textbf{29} (2) (1922), 51-55.

\bibitem{mr} M.~Rosiu;
 Trajectory structure near critical points,
 \emph{An. Univ. Craiova Ser. Mat. Inform.} \textbf{25} (1998), 35-44.

\bibitem{sps} G.~Stefanidou, G.~Papaschinopoulos, C.~Schinas;
On a system of max difference equations, \emph{Dynam. Contin.
Discrete Impuls. Systems Ser. A}, \textbf{14} (6) (2007), 885-903.

\bibitem{ps-cana} G. Stefanidou, G. Papaschinopoulos, C. J. Schinas;
 On a system of two exponential type difference equations,
 \emph{Commun. Appl. Nonlinear Anal.} \textbf{17} (2) (2010), 1-13.

\bibitem{ejde-fodce} S.~Stevi\' c;
 First-order product-type systems of difference equations solvable in closed form;
 \emph{Electron. J. Differential Equations} Vol. 2015, Article No. 308, 
(2015), 14 pages.

\bibitem{ejqtde-newm} S.~Stevi\' c;
 New solvable class of product-type systems of difference equations on the
 complex domain and a new method for proving the solvability,
 \emph{Electron. J. Qual. Theory Differ. Equ.} Vol. 2016, Article No. 120, 
(2016), 19 pages.

\bibitem{ana2} S.~Stevi\' c;
 Solvability of a product-type system of difference equations with six parameters,
 \emph{Adv. Nonlinear Anal.} (in press), doi:10.1515/anona-2016-0145.

\bibitem{ana1} S.~Stevi\' c;
 Solvable subclasses of a class of nonlinear second-order difference equations,
 \emph{Adv. Nonlinear Anal.} \textbf{5} (2) (2016), 147-165.

\bibitem{ejqtde-4th-2017} S.~Stevi\' c;
 Solvable product-type system of difference equations
whose associated polynomial is of the fourth order,
 \emph{Electron. J. Qual. Theory Differ. Equ.} Vol. 2017, Article No. 13, 
(2017), 29 pages.

\bibitem{ejde-2015} S.~Stevi\' c, M.~A.~Alghamdi, A.~Alotaibi, E.~M.~Elsayed;
 Solvable product-type system of difference equations of second order,
 \emph{Electron. J. Differential Equations} Vol. 2015, Article No. 169, (2015), 
 20 pages.

\bibitem{jia2015-pts} S.~Stevi\' c, B.~Iri\v{c}anin, Z. \v{S}marda;
On a product-type system of difference equations of second order solvable in 
closed form, \emph{J. Inequal. Appl.} Vol. 2015, Article No. 327, (2015), 15 pages.

\bibitem{ejde-2016} S.~Stevi\' c, B.~Iri\v{c}anin, Z.~\v{S}marda;
 Solvability of a close to symmetric system of difference equations; 
\emph{Electron. J. Differential Equations} Vol. 2016, Article No. 159, (2016),
 13 pages.

\bibitem{ade-16-ptdce} S.~Stevi\' c, B.~Iri\v{c}anin, Z.~\v{S}marda;
 Two-dimensional product-type system of difference equations solvable in
 closed form; \emph{Adv. Difference Equ.} Vol. 2016, Article No. 253, 
(2016), 20 pages.

\bibitem{ejqtde-st2016} S.~Stevi\' c and D.~Rankovi\'c;
 On a practically solvable product-type system of difference equations of 
second order, \emph{Electron. J. Qual. Theory Differ. Equ.} Vol. 2016, 
Article No. 56, (2016), 23 pages.

\end{thebibliography}

\end{document}