\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2013 (2013), No. 252, pp. 1--12.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2013/252\hfil Bogdanov-Takens bifurcation for NFDEs] {Bogdanov-Takens bifurcation for neutral functional differential equations} \author[J. Cao, R. Yuan \hfil EJDE-2013/252\hfilneg] {Jianzhi Cao, Rong Yuan} % in alphabetical order \address{Jianzhi Cao \newline School of Mathematical Sciences, Beijing Normal University Beijing 100875, China} \email{cjz2004987@163.com} \address{Rong Yuan \newline School of Mathematical Sciences, Beijing Normal University Beijing 100875, China} \email{ryuan@bnu.edu.cn} \thanks{Submitted August 6, 2013. Published November 20, 2013.} \thanks{Supported by the NSFC and RFDP} \subjclass[2000]{34K06, 34K18, 34K20, 34K60, 37G05, 37G10} \keywords{Neutral functional differential equations; center manifold; \hfill\break\indent Bogdanov-Takens bifurcation; normal forms} \begin{abstract} In this article, we consider a class of neutral functional differential equations (NFDEs). First, some feasible assumptions and algorithms are given for the determination of Bogdanov-Takens (B-T) singularity. Then, by employing the method based on center manifold reduction and normal form theory due to Faria and Magalh\~aes \cite{faria2}, a concrete reduced form for the parameterized NFDEs is obtained and the bifurcation behavior of the parameterized NFDEs is described. This result extend the B-T bifurcation analysis reported in \cite{xu}. Finally, two examples ilustrate the theoretical results. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{example}[theorem]{Example} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} As we all know, neutral functional differential equations (NFDEs) can be found in many scientific and technological fields such as biology, physics, control theory, engineering and so on. In the last decade, stability theory and Hopf bifurcation of NFDEs have been widely discussed by many authors because of its important(see \cite{babram, hale, han, niu2, qu, wang1,wang2, weeder} and the references therein). For example, in \cite{han}, by employing a Lyapunov functional, Han studies the following linear neutral system $\dot{x}(t) - C\dot{x}(t - \tau) = A(t)x(t) + B(t)x(t-r)$ with initial condition $x(t) = \varphi(t), \dot{x}(t) = \dot{\varphi} (t)$; in \cite{wang2}, the authors investigate a neutral delay logistic differential equation $\frac{d}{dt}[x(t)-cx(t - \tau)] = r(1 - e^{x(t-\sigma)})$ and derive explicit formulas for determining the direction of Hopf bifurcation; Qu et al~\cite{qu} prove that the following equation $\frac{d}{dt}[x(t)-px(t - \tau)] = -ax(t)+b\tanh( x(t-\tau ))$ has a sequence of Hopf bifurcations. Moreover, they illustrate the global existence of periodic solutions by using global Hopf bifurcation technique; In \cite{weeder}, using center manifold reduction and normal form theory, Weedermann derive explicit conditions for a Hopf bifurcation in scalar neutral system $\dot{x}(t)-a\dot{x}(t-\tau)=L(\gamma)x_t+f(x_t,\gamma).$ Recently, codimension 2 bifurcations, such as B-T bifurcation, zero-Hopf bifurcation and so on, become a subject of intense research activities (see \cite{jiang, niu1, niu2, xu, yuan}). The phenomenon of B-T bifurcation is a bifurcation of an equilibrium point in a two-parameter family of autonomous differential equations at which the critical equilibrium has a zero eigenvalue of (algebraic) multiplicity two. For nearby parameter values, the system has two equilibria (a saddle and a nonsaddle) which collide and disappear via a saddle-node bifurcation. The nonsaddle equilibrium undergoes an Hopf bifurcation generating a limit cycle. This cycle degenerates into an orbit homoclinic to the saddle and disappears via a saddle homoclinic bifurcation (see Figure \ref{fig1}). As far as we know, till now, few works are concerned with the B-T bifurcations of the NFDEs(see \cite{niu1}). This fact motivates our work for the manuscript. The paper is organized as follows: in section 2, we characterize the B-T singularity for parameterized NFDEs of the form \eqref{e2.1} in $\mathbb{R}^n$ with $n>1$, and give feasible algorithms to calculate the explicit expression of the generalized eigenspace associated with the zero eigenvalue of the linearized NFDEs; in section 3, by employing the methods used in \cite{faria1, faria2}, we reduce the parameterized NFDEs to the normal form on a center manifold, and a concrete expression of the reduced form is obtained. On the other hand, the bifurcation properties caused by B-T singularity are described in Proposition \ref{prop3.2}. In the last section, two examples are presented to illustrate our results. \section{Theory and algorithm for NFDEs} First, we consider the system $$\label{e2.1} \frac{d}{dt}\mathcal {D}u_t=\mathcal {L}(\alpha)u_t+\mathcal {F}(\alpha,u_t)$$ where $u_t\in C_n:=C([-r, 0],\mathbb{R}^n)$, $\alpha\in\mathbb{R}^2$, $\mathcal {D}$ and $\mathcal{L}(\alpha)$ are continuous linear functions from $C_n$ to $\mathbb{R}^n$, $\mathcal {F}(\alpha,0)=\frac{\partial\mathcal {F}}{\partial x}(\alpha,0)=0$. Suppose that $\mu$ and $\eta$ are $n\times n$ matrix valued functions of bounded variation in $\theta\in[-r,0]$ and $\mu$ is non-atomic at zero (see \cite{hale} for the definition of atomic at zero) and define \begin{gather*} \mathcal {D}\phi=\phi(0)-\int_{-r}^0d[\mu(\theta)]\phi(\theta),\\ \mathcal {L}(\alpha)\phi=\int_{-r}^0d[\eta(\alpha,\theta)]\phi(\theta), \end{gather*} for all $\phi\in C_n$ . We consider the linearized equation of \eqref{e2.1}, $$\label{e2.2} \frac{d}{dt}\mathcal {D}u_t=\mathcal {L}(\alpha)u_t.$$ If we denote by $x(.,\phi)$ the unique solution of \eqref{e2.2} with initial function $\phi$ at zero, then \eqref{e2.2} determines a $C_0-$semigroup of bounded linear operators given by $T(t)\phi=x_t(\phi),\quad \text{for } t\geq 0,$ where $x_t(\phi)$ is the solution of \eqref{e2.2} with $x_0(\phi)=\phi$. Denote by $\mathcal {A}$ the infinitesimal generator of $\{T(t)\}_{t\geq 0}$ $$\label{e2.3} D(\mathcal {A})=\Big\{\phi\in C_n: \frac{d\phi}{d\theta} \in C_n, \mathcal {D}\frac{d\phi}{d\theta}=\mathcal {L}(\alpha)\phi\Big\},\quad \mathcal {A}\phi=\frac{d\phi}{d\theta}.$$ If $\mathcal {A}$ is defined by \eqref{e2.3}, then $\sigma(\mathcal {A})=P\sigma(\mathcal {A})$ and $\lambda$ is in $\sigma(\mathcal {A})$ if and only if satisfies the characteristic equation $$\label{e2.4} \det \Delta(\lambda)=0,\ \ \Delta(\lambda) =\lambda\mathcal {D}(e^{\lambda.}I)-\mathcal {L}(\alpha)(e^{\lambda.}I).$$ Further assumptions on system \eqref{e2.1} are as follows: \begin{itemize} \item[(H1)] $\lambda=0$ is a characteristic value of $\mathcal {A}_0$ with algebraic multiplicity 2 and geometric multiplicity 1, where $\mathcal {A}_0$ is the generator of $\{T(t)\}_{t\geq 0}$ when $\alpha=0$; \item[(H2)] all other characteristic values of $\mathcal {A}_0$ have non-zero real parts. \end{itemize} We say that system \eqref{e2.1} has a B-T singularity if (H1)-(H2) hold, and in this case we call $(x,\alpha)=(0,0)$ a B-T point of system \eqref{e2.1}. In this article, we consider the system $$\label{e2.5} \frac{d}{dt}(x(t)+Ex(t-r))=A(\alpha)x(t)+B(\alpha)x(t-r)+F(\alpha,x(t),x(t-r))$$ where $E, A(\alpha), B(\alpha)\in \mathbb{R}^{n\times n}$, $x\in \mathbb{R}^n$, and $F(\alpha,0,0)=\frac{\partial F(\alpha,0,0)}{\partial x} =\frac{\partial F(\alpha,0,0)}{\partial y}=0$. If we choose $\mu(\theta)= \begin{cases} 0 & -r<\theta\leq 0, \\ E & \theta=-r, \end{cases}$ and $\eta(\alpha,\theta)= \begin{cases} A(\alpha)+B(\alpha) & \theta=0, \\ B(\alpha) & -r<\theta<0, \\ 0 & \theta=-r, \end{cases}$ then system \eqref{e2.5} can be formulated in the form \eqref{e2.1}; i.e. $$\label{e2.6} \frac{d}{dt}\mathcal {D}x_t=\mathcal {L}(\alpha)x_t+\mathcal {F}(\alpha,x_t).$$ Further, we have $$\label{e2.7} \frac{d}{dt}\mathcal {D}x_t=\mathcal {L}_0x_t+\tilde{\mathcal {F}}(\alpha,x_t),$$ where $\mathcal {L}_0:=\mathcal {L}(0)$ and $\tilde{\mathcal {F}}(\alpha,x_t)=(\mathcal {L}(\alpha) -\mathcal {L}_0)x_t+\mathcal {F}(\alpha,x_t)$. Obviously, the linearized equation of system \eqref{e2.7} at $(0,0)$ is $$\label{e2.8} \frac{d}{dt}\mathcal {D}x_t=\mathcal {L}_0x_t.$$ The following theorem gives an equivalent description for B-T singularity in a neutral differential system \eqref{e2.5}, which can be used as a feasible algorithm for determining the B-T singularity. \begin{theorem}\label{thm2.1} Under assumption (H2), system \eqref{e2.5} has a B-T singularity if and only if the following conditions hold: \begin{itemize} \item[(i)] $\operatorname{rank}(A+B)=n-1$; \item[(ii)] if $\ker (A+B)=\operatorname{span}\{\phi_1^0\}$, then $(I+E+rB)\phi_1^0\in\operatorname{Ran}(A+B)$; \item[(iii)] if $(I+E+rB)\phi_1^0=(A+B)\phi_2^0$, then $(I+E+rB)\phi_2^0-(rE+\frac{1}{2}r^2B)\phi_1^0 \notin\operatorname{Ran}(A+B)$, where $\phi_1^0, \phi_2^0\in \mathbb{R}^n$, $A=A(0)$ and $B=B(0)$. \end{itemize} \end{theorem} \begin{proof} Similar to the reference \cite{xu}, we only needs to show the equivalence of assumption (H1) and the conditions (i)--(iii). Here we translate assumption (H1) as follows: there exist linearly independent functions $\phi_1,\ \phi_2\in C_n$ such that $$\label{e2.9} \mathcal {A}_0\phi_1=0,\ \ \mathcal {A}_0\phi_2=\phi_1$$ and the following equation $$\label{e2.10} \mathcal {A}_0\phi=\phi_2$$ has no solution $\phi$ in $C_n$. Notice that $\mathcal {A}_0\phi_1=0$ is equivalent to \begin{gather*} \dot{ \phi}_1(\theta)=0 -r\leq \theta< 0,\\ \mathcal {D}\dot{\phi}_1=\mathcal {L}_0\phi_1 \quad\theta=0, \end{gather*} from which we obtain $$\label{e2.11} \phi_1(\theta)\equiv \phi_1^0\in \mathbb{R}^n\setminus\{0\}\quad \text{and}\quad (A+B)\phi_1^0=0.$$ Moreover, from \eqref{e2.9} we see that the associated eigenspace of $\mathcal {A}_0$ is a space with dimension $2$ spanned by the characteristic function $\phi_1$ and the generalized characteristic function $\phi_2$, and that any function $\tilde{\phi}_1\in C_n$ linearly independent with $\phi_1$ must satisfy $\mathcal {A}_0\tilde{\phi}_1\neq 0,$ or $$\label{e2.12} (A+B)\tilde{\phi}_1^0\neq 0$$ for any $\tilde{\phi}_1^0$ in $\mathbb{R}^n$ linearly independent with $\phi_1^0$. As a consequence of \eqref{e2.11} and \eqref{e2.12}, the condition (i) holds; i.e., $\operatorname{rank}(A+B)=n-1$ and $\ker(A+B)=\operatorname{span}\{\phi_1^0\}$. By the definition of the operator $\mathcal {A}_0$ and \eqref{e2.11}, $\mathcal {A}_0\phi_2=\phi_1$ is equivalent to $$\label{e2.13} \begin{gathered} \dot{\phi}_2(\theta)=\phi_1^0,\quad -r\leq\theta<0,\\ \mathcal {D}\dot{\phi}_2=\mathcal {L}_0\phi_2,\quad \theta=0, \end{gathered}$$ which implies $$\label{e2.14} \phi_2(\theta)=\phi_2^0+\phi_1^0\theta,$$ and $\phi_1^0+E\phi_1^0=(A+B)\phi_2^0-rB\phi_1^0;$ i.e., $(A+B)\phi_2^0=(I+E+rB)\phi_1^0$, the condition (ii) holds. Similar to the above discussion, \eqref{e2.10} is equivalent to $$\label{e2.15} \begin{gathered} \dot{\phi}(\theta)=\phi_2^0+\phi_1^0\theta,\quad -r\leq\theta<0,\\ \mathcal {D}\dot{\phi}=\mathcal {L}_0\phi,\quad \theta=0. \end{gathered}$$ From the first formula of \eqref{e2.15}, we have $\phi(\theta)=\phi^0+\phi_2^0(\theta)+\frac{1}{2}\phi_1^0\theta^2$, $\phi^0\in \mathbb{R}^n$, and by substituting $\phi(\theta)$ into the second formula, we obtain $$\label{e2.16} (I+E+rB)\phi_2^0-(rE+\frac{1}{2}r^2B)\phi_1^0=(A+B)\phi^0.$$ From \eqref{e2.16} we see that \eqref{e2.10} has no solution in $C_n$ is equivalent to $(I+E+rB)\phi_2^0-(rE+\frac{1}{2}r^2B)\phi_1^0\notin\operatorname{Ran}(A+B)$; i.e. the condition (iii) holds. This completes the proof of the theorem. \end{proof} We can easily obtain the following corollaries. \begin{corollary}\label{cor2.2} The conditions {\rm (ii)} and {\rm (iii)} in Theorem \ref{thm2.1} are equivalent to \begin{itemize} \item[(ii')] if $\ker (A+B)=\operatorname{span}\{\phi_1^0\}$, then $\operatorname{rank}(A+B,(I+E+rB)\phi_1^0)=n-1$; \item[(iii')] if $(I+E+rB)\phi_1^0=(A+B)\phi_2^0$, then $\operatorname{rank}(A+B,(I+E+rB)\phi_2^0-(rE+\frac{1}{2}r^2B)\phi_1^0)=n$. \end{itemize} \end{corollary} \begin{corollary}\label{cor2.3} When $n=1$ (i.e. $E, A, B\in \mathbb{R}$), Equation \eqref{e2.5} has a B-T singularity if and only if $A=-B$, $1+E=-rB$ and $E\neq1$. \end{corollary} In the space $C_n^*=C([0,r],\ \mathbb{R}^{n*})$, the adjoint bilinear form on $C_n^*\times C_n$ is defined by $(\psi,\phi)=\psi(0)\phi(0)-\int_{-r}^0d[\int_0^\theta\psi(\theta-\xi) d[\mu(\xi)]]\phi(\theta) -\int_{-r}^0\int_0^\theta\psi(\theta-\tau)d[\eta(\alpha,\tau)]\phi(\theta).$ For the set $\Lambda=\{0\}$, we can use the formal adjoint theory for an NFDE to decompose the phase space. The hypotheses (H1)--(H2) yield a decomposition of the state space in the form $C_n=P\oplus Q$, where $P$ is the two dimensional generalized eigenspace associated to zero for \eqref{e2.8} and $Q=\{\phi\in C_n:(\psi,\phi)=0$\} is the complementary space. Let $\Phi(\theta)=(\phi_1(\theta),\phi_2(\theta))$, $-r\leq\theta\leq0$ and $\Psi(s)=\operatorname{col}(\psi_1(s),\psi_2(s))$, $0\leq s\leq r$ be the bases of $P$ and its dual space $P^*$, respectively, which satisfy $(\Psi,\Phi)=I_2$($I_2$ is the $2\times 2$ identity matrix). Then, we have the following lemma. \begin{lemma}\label{lem2.4} The bases of $P$ and its dual space $P^*$ have the following representations \begin{gather} P=\operatorname{span}\Phi, \quad \Phi(\theta)=(\phi_1(\theta),\phi_2(\theta)), \quad -r\leq\theta\leq0, \\ P^*=\operatorname{span}\Psi \quad \Psi(s)=col(\psi_1(s),\psi_2(s)), \quad 0\leq s\leq r , \end{gather} where $\phi_1(\theta)= \phi_1^0\in \mathbb{R}^n\setminus\{0\}$, $\phi_2(\theta)=\phi_2^0+\phi_1^0\theta$, $\phi_2^0\in \mathbb{R}^n$, and $\psi_2(s)=\psi_2^0\in\mathbb{R}^{n*}\setminus\{0\}$, $\psi_1(s)=\psi_1^0-s\psi_2^0$, $\psi_1^0\in\mathbb{R}^{n*}$, which satisfy \begin{itemize} \item[(i)] $(A+B)\phi_1^0=0$, \item[(ii)] $(A+B)\phi_2^0=(I+E+rB)\phi_1^0$, \item[(iii)] $\psi_2^0(A+B)=0$, \item[(iv)] $\psi_1^0(A+B)=\psi_2^0(I+E+rB)$, \item[(v)] $\psi_2^0(I+E+rB)\phi_2^0-\psi_2^0(rE+\frac{1}{2}r^2)B\phi_1^0=1$, \item[(vi)] $\psi_1^0(I+E+rB)\phi_2^0-\psi_1^0(rE+\frac{1}{2}r^2B)\phi_1^0 +\psi_2^0(r^2E+\frac{1}{6}r^3B)\phi_1^0-\psi_2^0(rE+\frac{1}{2}r^2B)\phi_2^0=0$. \end{itemize} Here, we can fix the unique vectors $\phi_1^0,\ \phi_2^0$ by (i), (ii), respectively; then $\psi_2^0,\ \psi_1^0$ can be determined by (iii), (iv), respectively, up to some constant factors. Furthermore, (v) and (vi) are used to determine the coefficient factors of the vectors $\psi_2^0$ and $\psi_1^0$. \end{lemma} \begin{proof} Firstly, from the proof of Theorem \ref{thm2.1}, we know that $\phi_1(\theta)= \phi_1^0\in \mathbb{R}^n\setminus\{0\}$, $\phi_2(\theta)=\phi_2^0+\phi_1^0\theta$, $\phi_2^0\in \mathbb{R}^n$, and then (i), (ii) hold. Next, we introduce the adjoint operator $\mathcal {A}_0^*: C_n^*\to C_n^*$ of $\mathcal {A}_0$ by \begin{gather*} D(\mathcal {A}_0^*)=\Big\{ \psi\in C^1([0,r],\mathbb{R}^{n*}) =:C_n^{*1}:\frac{d\psi}{ds}\in C_n^{*1},\; \mathcal {D}\frac{d\psi}{ds} =- \int_{-r}^0\psi(-\theta)\textrm{d}\eta(\theta)\Big\}, \\ \mathcal {A}_0^*=-\frac{d\psi}{ds}. \end{gather*} From $\mathcal {A}_0^*\psi_2=0$, i.e. $$\label{e2.17} \begin{gathered} -\dot{\psi}_2(s)=0,\quad 0< s\leq r,\\ \mathcal {D}\dot{\psi}_2=-\mathcal {L}_0\psi_2,\quad s=0,\\ \end{gathered}$$ we obtain $$\label{e2.18} \psi_2(s)=\psi_2^0\in\mathbb{R}^{n*}\setminus\{0\},\quad \psi_2^0(A+B)=0.$$ From $\mathcal {A}_0^*\psi_1=\psi_2$, i.e. $$\label{e2.19} \begin{gathered} -\dot{\psi}_1(s)=\psi_2^0,\quad 0< s\leq r,\\ \mathcal {D}\dot{\psi}_1=-\mathcal {L}_0\psi_1,\quad s=0, \end{gathered}$$ we have $$\label{e2.20} \psi_1(s)=\psi_1^0-s\psi_2^0, \quad \psi_1^0\in\mathbb{R}^n,\quad \psi_1^0(A+B)=\psi_2^0(I+E+rB)$$ which implies that the conditions (iii) and (iv) hold. Finally, by the definition of $(\cdot,\cdot)$ and ($\Phi,\Psi)=1$ we have \label{e2.21} \begin{gathered} (\psi_1,\phi_1)=\psi_1^0(I+E+rB)\phi_1^0-\psi_2^0(rE+\frac{1}{2}r^2)B\phi_1^0=1, \\ (\psi_2,\phi_2)=\psi_2^0(I+E+rB)\phi_2^0-\psi_2^0(rE+\frac{1}{2}r^2)B\phi_1^0=1, \\ \begin{aligned} (\psi_1,\phi_2)&=\psi_1^0(I+E+rB)\phi_2^0-\psi_1^0(rE+\frac{1}{2}r^2B)\phi_1^0\\ &\quad +\psi_2^0(r^2E+\frac{1}{6}r^3B)\phi_1^0 -\psi_2^0(rE+\frac{1}{2}r^2B)\psi_2^0=0, \end{aligned} \\ (\psi_2,\phi_1)=\psi_2^0(I+E+rB)\phi_1^0=0. \end{gathered} In fact, according to the above formulas, we know that the fourth formula in \eqref{e2.21} holds naturally, and the first is equivalent to the second one. Then we can properly choose the coefficient factors of $\psi_1^0$, $\psi_2^0$ such that all the formulas hold. This completes the proof of the lemma. \end{proof} \section{Reduction and normal form for NFDEs with B-T bifurcation} In this section, we will discuss the reduction and normal forms for the NFDEs with B-T singularity. As we all know, one way of considering normal forms for NFDEs with parameters is to reduce the situation to the case of differential equations without parameters by considering the system $$\label{e3.1} \begin{gathered} \frac{d}{dt}\mathcal {D}x_t=\mathcal {L}_0x_t+\tilde{\mathcal {F}}(\alpha,x_t),\\ \frac{d}{dt}\alpha(t) =0. \end{gathered}$$ The solutions of system \eqref{e3.1} are of the form $\tilde{x}(t)=(x(t),\alpha(t))\in\mathbb{R}^n\times\mathbb{R}^2$ and the phase space is $C_{n+2}:=C([-r,0],\mathbb{R}^n\times\mathbb{R}^2)$. Let $\tilde{\mathcal {L}}_0\tilde{x}_t=(\mathcal {L}_0x_t,0)$ and $\hat{\mathcal {F}}\big(\tilde{x}_t)=(\tilde{\mathcal {F}}(\alpha(0),x_t),0)$. It is now possible to apply to \eqref{e3.1} the normal form theory for FDEs without parameters developed in \cite{faria1}. Denote $\tilde{\mathcal {A}}_0=(\mathcal {A}_0,0)$ then the eigenvalues of $\tilde{\mathcal {A}}_0$ include not only all the eigenvalues of $\mathcal {A}_0$, but also the two $0$-eigenvalues introduced by $\dot{\alpha}=0$. Let $\tilde{\Lambda}$ be the set of $0$-eigenvalues (counting multiplicity) of $\tilde{\mathcal {A}}_0$. We consider the decomposition of $C_{n+2}=\tilde{P}\oplus \tilde{Q}$, where $\tilde{P}=P\times\mathbb{R}^2$ is the invariant space of $\tilde{\mathcal {A}}_0$ associated with $\tilde{\Lambda}$, $\tilde{Q}=Q\times R$, where $R=\{v\in C_2\mid v(0)=0\}$. Here, the bases of $\tilde{P}$ and its dual space $\tilde{P}^*$ are formed by the columns of $\tilde{\Phi}= \begin{bmatrix} \Phi& 0 \\ 0 & I_2 \end{bmatrix}$ and the rows of $\tilde{\Psi}= \begin{bmatrix} \Psi & 0 \\ 0 & I_2 \end{bmatrix}$, respectively, which satisfy $(\tilde{\Psi},\tilde{\Phi})=I_4$ and $\dot{\tilde{\Phi}}=\tilde{\Phi}\tilde{J}$ with $\tilde{J}= \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}.$ Let $X_0$ and $Y_0$ be the matrix valued functions defined on $[-r,0]$ by \begin{gather*} X_0(\theta)= \begin{cases} 0 & -r\leq\theta<0 \\ I_n &\theta=0 , \end{cases} \\ Y_0(\theta)= \begin{cases} 0 & -r\leq\theta<0 \\ I_2 &\theta=0\,. \end{cases} \end{gather*} As shown in \cite{faria1}, we consider the enlarged phase space $BC_n$ in which the functions from $[-r,0]$ to $\mathbb{R}^n$ are uniformly continuous on $[-r,0)$ with a jump discontinuously at $0$, this space can be identified as $C_n\times \mathbb{R}^n$; its elements are in the form $\phi=\varphi+X_0\xi$, where $\varphi\in C_n$, $\xi\in\mathbb{R}^n$. Analogously, for considering normal forms for the equation \eqref{e3.1} we take for phase space $B\tilde{C}_n:=BC_n\times BC_2$, which can be identified with $C_{n+2}\times\mathbb{R}^{n+2}$. The definition of the continuous projection $\pi:\ BC_n\to P, \quad \pi(\phi+X_0\xi)=\Phi[(\Psi,\phi)+\Psi(0)\xi]$ allows us to decompose the enlarged phase space by $BC_n=P\oplus$ $\ker \pi$ and $Q\subset\ker \pi$. In the space $BC_n$, \eqref{e2.5} can be written as an abstract ordinary differential equation $$\label{e3.2} \frac{d}{dt}x_t=\mathcal {A}_0x_t+X_0\tilde{\mathcal{F}}(\alpha,x_t),$$ where $\mathcal{A}_0\phi=\phi'+X_0[\mathcal{L}_0\phi-\mathcal{D}\phi']$. We consider the projection $\tilde{\pi}: B\tilde{C}_n\to\tilde{P}$ given by $\tilde{\pi}(\phi+X_0\xi,\psi+Y_0\zeta)=\tilde{\Phi} \Big\{\Big(\tilde{\Psi}, \begin{bmatrix} \phi \\ \psi \end{bmatrix} \Big)+\tilde{\Psi}(0) \begin{bmatrix} \xi \\ \zeta \end{bmatrix}\Big\} =(\pi(\phi+X_0\xi),\psi(0)+\zeta)$ where $\phi\in C_n,\ \psi\in C_2,\ \xi\in\mathbb{R}^n$ and $\zeta\in\mathbb{R}^2$. According to $B\tilde{C}_n=\tilde{P}\oplus \ker \tilde{\pi}$, we can decompose $\begin{bmatrix} x_t \\ \alpha_t \end{bmatrix} = \begin{bmatrix} \Phi & 0 \\ 0 & I_2 \end{bmatrix} \begin{bmatrix} x(t) \\ \alpha(t) \end{bmatrix} + \begin{bmatrix} y \\ w \end{bmatrix},$ where $(x(t),\alpha(t))\in \mathbb{R}^{n+2}$ and $(y,w)\in \ker \tilde{\pi}$. Thus $$\label{e3.3} \begin{gathered} \begin{bmatrix} \dot{x} \\ \dot{\alpha} \end{bmatrix} =\tilde{J} \begin{bmatrix} x \\ \alpha \end{bmatrix} +\tilde{\Psi}(0)\hat{\mathcal {F}}\Big( \tilde{\Phi} \begin{bmatrix} x \\ \alpha \end{bmatrix} + \begin{bmatrix} y \\ w \end{bmatrix} \Big) \\ \frac{d}{dt} \begin{bmatrix} y \\ w \end{bmatrix} = \tilde{\mathcal {A}}_{Q^1} \begin{bmatrix} y \\ w \end{bmatrix} +(I-\pi)[X_0,Y_0]\hat{\mathcal {F}} \Big( \tilde{\Phi} \begin{bmatrix} x \\ \alpha \\ \end{bmatrix} + \begin{bmatrix} y \\ w \end{bmatrix} \Big) \end{gathered}$$ where $x\in \mathbb{R}^2$, $\alpha\in\mathbb{R}^2$, $y\in Q^1:=Q\cap C_n^1$, $w\in R^1:=R\cap C_2^1$, and $\tilde{\mathcal {A}}_{Q^1}$ is the operator from $\tilde{Q}^1:=\tilde{Q}\cap C_{n+2}^1=Q^1\times R^1$ to $\ker\tilde{\pi}$, $\tilde{\mathcal {A}}_{\tilde{Q}}^1 \begin{bmatrix} \phi \\ \psi \end{bmatrix} = \begin{bmatrix} \dot{\phi} \\ \dot{ \psi} \end{bmatrix} +[X_0,Y_0] \Big\{ \tilde{\mathcal {L}}_0 \begin{bmatrix} \phi \\ \psi \end{bmatrix} -\mathcal {D} \begin{bmatrix} \dot{\phi}(0) \\ \dot{\psi}(0) \end{bmatrix} \Big\}.$ If $\mathcal{A}_{Q^1}\subset\ker \pi$ is such that $\mathcal {A}_{Q}^1\phi=\dot{\phi}+X_0[\mathcal {L}_0\phi -\mathcal {D}\dot{\phi}(0)]$, system \eqref{e3.3} is equivalent to \begin{align*} \begin{bmatrix} \dot{x} \\ \dot{\alpha} \end{bmatrix} &= \begin{bmatrix} Jx \\ 0 \end{bmatrix} + \begin{bmatrix} \Psi(0)(\mathcal {L}(\alpha(0))+w(0)-\mathcal {L}_0)(\Phi x+y) \\ 0 \end{bmatrix} \\ &\quad + \begin{bmatrix} \Psi(0)\mathcal {F}(\Phi x+y,\alpha(0)+w(0)) \\ 0 \end{bmatrix}, \end{align*} \begin{align*} \begin{bmatrix} \dot{y} \\ \dot{w} \end{bmatrix} &= \begin{bmatrix} \mathcal {A}_{Q^1}y \\ \dot{w}-Y_0\dot{w}(0) \end{bmatrix} +\begin{bmatrix} (I-\pi)X_0(\mathcal {L}(\alpha(0))+w(0))-\mathcal {L}_0)(\Phi x+y) \\ 0 \end{bmatrix} \\ &\quad + \begin{bmatrix} (I-\pi)X_0\mathcal{F}(\Phi x+y,\alpha(0)+w(0)) \\ 0 \end{bmatrix}. \end{align*} Noting that $w(0)=0$, system \eqref{e3.3} can be reduced to the following equation on $BC_n$, $$\label{e3.4} \begin{gathered} \dot{x}=Jx+\Psi(0)\tilde{\mathcal{F}}(\Phi x+y,\alpha) \\ \frac{d}{dt}y= \mathcal {A}_{Q^1}y+(I-\pi)X_0\tilde{\mathcal{F}} (\Phi x+y,\alpha). \end{gathered}$$ Employing Taylor's theorem, we expand the nonlinear terms in \eqref{e3.4} at $(x,y,\alpha)=(0,0,0)$ as $$\label{e3.5} \begin{gathered} \dot{x}=Jx+\sum_{j\geq2}\frac{1}{j!}f_j^1(x,y,\alpha), \\ \frac{d}{dt} y=\mathcal {A}_{Q^1}y+\sum_{j\geq2}\frac{1}{j!}f_j^2(x,y,\alpha), \end{gathered}$$ where $f_j^i(x,y,\alpha)(i=1,2)$ denotes the homogeneous polynomials of degree $j$ in variables $(x,y,\alpha)$. Let $V_2^4(\mathbb{R}^2)= \big\{ \sum_{|(q,l)=2|}c _{(q,l)}x^q\alpha^l:(q,l)\in\mathbb{N}_0^4,\; c_{(q,l)}\in\mathbb{R}^2 \big\}$ denote the linear space of homogeneous polynomials of $(x,\alpha)=(x_1,x_2,\alpha_1,\alpha_2)$ with the degree 2 and coefficients in $\mathbb{R}^2$ and let $M_2^1$ be the operator defined on $V_2^4(\mathbb{R}^2)$ by $(M_2^1p)(x,\alpha)=D_xp(x,\alpha)Jx-Jp(x,\alpha),\quad \forall p\in V_2^4(\mathbb{R}^2),$ then we decompose $V_2^4(\mathbb{R}^2)$ as $\operatorname{Im}(M_2^1)\oplus \operatorname{Im}(M_2^1)^c$ and denote the map from $V_2^4(\mathbb{R}^2)$ to $\operatorname{Im}(M_2^1)$ by $P_{I,2}^1$. We start from finding the second-order normal form($j=2$). The canonical basis of $V_2^4(\mathbb{R}^2)$ has twenty elements $((x,\alpha)^2,0)^T,\quad (0,(x,\alpha)^2)^T$ and we can choose the complementary space $\operatorname{Im}(M_2^1)^c$ spanned by the elements $\begin{bmatrix} 0 \\ x_1^2 \end{bmatrix},\quad \begin{bmatrix} 0 \\ x_1x_2 \end{bmatrix},\quad \begin{bmatrix} 0 \\ x_1\alpha_i \end{bmatrix},\quad \begin{bmatrix} 0 \\ x_2\alpha_i \end{bmatrix}, \quad \begin{bmatrix} 0 \\ \alpha_i^2 \end{bmatrix},\quad \begin{bmatrix} 0 \\ \alpha_1\alpha_2 \end{bmatrix},\quad i=1,2.$ Then, for the system \eqref{e2.5}, the normal form with universal unfolding of double zero singularity has the form $\dot{x}=Jx+\frac{1}{2}g_2^1(x,0,\alpha)+\text{h.o.t.,}$ where $g_2^1(x,0,\alpha)=(I-P_{I,2}^1)f_2^1(x,0,\alpha) =\operatorname{proj}_{Im(M_2^1)^c}f_2^1(x,0,\alpha)$ is the quadratic function in $(x,\alpha)$ for $y=0$, and h.o.t. stands for higher order terms. Following \cite{faria2, xu}, the normal form of system \eqref{e3.1} on the center manifold is given by $$\label{e3.6} \dot{x}_1 =x_2,\quad \dot{x}_2=\rho_1x_1+\rho_2x_2+\eta_1x_1^2+\eta_2x_1x_2+\text{h.o.t.,}$$ where \begin{gather*} \rho_1=\psi_2^0(A_1+B_1)\phi_1^0\alpha_1+\psi_2^0(A_2+B_2)\phi_1^0\alpha_2, \\ \begin{aligned} \rho_2&=[\psi_1^0(A_1+B_1)\phi_1^0+\psi_2^0((A_1+B_1)\phi_2^0-B_1\phi_1^0)] \alpha_1 \\ &\quad +[\psi_1^0(A_2+B_2)\phi_1^0+\psi_2^0((A_2+B_2)\phi_2^0-B_2\phi_1^0)] \alpha_2, \end{aligned} \\ \eta_1=\psi_2^0\sum_{i=1}^n(E_i+F_i+G_i)\phi_1^0\phi_{1i}^0 ,\\ \begin{aligned} \eta_2&=2 \psi_1^0\sum_{i=1}^n(E_i+F_i+G_i)\phi_1^0\phi_{1i}^0\\ &\quad +\psi_2^0\Big[\sum_{i=1}^n(E_i+F_i+G_i)(\phi_2^0\phi_{1i}^0 +\phi_1^0\phi_{2i}^0) -\sum_{i=1}^n(E_i+2G_i)\phi_1^0\phi_{1i}^0\Big]. \end{aligned} \end{gather*} and $A_i=(A(\alpha)-A)'\big|_{\alpha_i}$, $B_i=(B(\alpha)-B)'\big|_{\alpha_i}$, $i=1,2$, $E_j, F_j, G_j$ $j=1,2,\ldots,n$ are coefficient matrices. \begin{remark} \rm Compared with \cite{xu}, we find that NFDEs and DDEs have the same normal forms in terms of B-T bifurcation except for the $\Phi(\theta)$ and $\Psi(s)$. \end{remark} After neglecting the terms of order higher than two in \eqref{e3.6}, we obtain $$\label{e3.7} \begin{gathered} \dot{x}_1 =x_2,\\ \dot{x}_2=\rho_1x_1+\rho_2x_2+\eta_1x_1^2+\eta_2x_1x_2. \end{gathered}$$ Next, we will list some dynamic behaviors about system \eqref{e3.7}. \begin{proposition}[{\cite[Theorem 4.3]{xu}}] \label{prop3.2} If $\eta_1\cdot \eta_2\neq0$ $(\eta_1, \eta_2\in\mathbb{R}$ defined in \eqref{e3.7}), then \eqref{e3.7} exhibits a generate B-T bifurcation. Furthermore, there exists a constant $\rho_1^0>0$, such that when $0<\rho_1(\alpha_1,\alpha_2)<\rho_1^0$, in the parameter plane ($\alpha_1,\alpha_2$) near the origin there exist two curve H and HL: \begin{itemize} \item The Hopf bifurcation curve $$H=\big\{(\alpha_1,\alpha_2):\rho_2(\alpha_1,\alpha_2) =\frac{\eta_2}{\eta_1}\rho_1(\alpha_1,\alpha_2)+\text{h.o.t.}=0,\; \rho_1(\alpha_1,\alpha_2)>0\big\};$$ \item The homoclinic bifurcation curve \begin{align*} HL=\Big\{&(\alpha_1,\alpha_2):\rho_2(\alpha_1,\alpha_2) -\mu(\sqrt{\rho_1(\alpha_1,\alpha_2)})\rho_1(\alpha_1,\alpha_2)+ \text{h.o.t.}=0,\\ &\rho_1(\alpha_1,\alpha_2)>0\Big\}; \end{align*} where $\mu$ is a continuously differentiable function with $\mu(0)=7\eta_1/(6\eta_2)$. \end{itemize} \end{proposition} To understand the propositions of intuitive, we give some numerical simulation (see Figure \ref{fig1} and Figure \ref{fig2}). \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{fig1a} % 12.eps \includegraphics[width=0.48\textwidth]{fig1b} \\ % 123.eps (a) \hfil (b) \end{center} \caption{For $\eta_1=-6/5$ and $\eta_2=2/5$, (a) shows that when $(\rho_1,\rho_2)=(0,0)$, the unique equilibrium $(0,0)$ is a cusp of codimension 2; (b) shows that when $(\rho_1,\rho_2)=(1.67,-0.934)$, there are two equilibria: one is a stable focus and the other is a saddle} \label{fig1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{fig2a} % 122.eps \includegraphics[width=0.48\textwidth]{fig2b} \\ % 132.eps (a) \hfil (b) \end{center} \caption{For $\eta_1=-6/5$ and $\eta_2=2/5$, (a) shows that when $(\rho_1,\rho_2)=(1.364,-0.446)$, the stable focus becomes unstable and there is a stable limit cycle, i.e., system \eqref{e3.6} undergoes a hopf bifurcation; (b) shows that when $(\rho_1,\rho_2)=(1.0833,1.6782)$, the limit cycle is broken and reach the manifold of the saddle, the homoclinic loop occur} \label{fig2} \end{figure} \section{Examples} For the applications of our method proposed in the previous section, we investigate B-T bifurcation exhibited in the following two examples. \begin{example}\label{exmp4.1}\rm We consider the following scalar neutral differential equation $$\label{e4.1} \dot{x}(t)+p\dot{x}(t-1)=a_1x(t)-a_2x(t-1)+ax^2(t)+bx(t)x(t-1)+cx^2(t-1)$$ where $a_1, a_2, a, b, c, p\in \mathbb{R}$. The characteristic equation for its linearization at zero is given by $\lambda(1+pe^{-\lambda})=a_1-a_2e^{-\lambda}.$ From Corollary \ref{cor2.3}, one can check that system \eqref{e4.1} has a B-T singularity if and only if $a_1=a_2,\quad 1+p=a_2,\quad p\neq 1.$ Rescaling the bifurcation parameters by setting $\alpha_1=a_1-1-p$, $\alpha_2=a_2-1-p$, $\alpha=(\alpha_1,\alpha_2)\in \mathbb{R}^2$, system \eqref{e4.1} becomes \begin{align*} \dot{x}(t)+p\dot{x}(t-1) &= (1+p)x(t)-(1+p)x(t-1)+\alpha_1x(t)-\alpha_2x(t-1) \\ &\quad +ax^2(t)+bx(t)x(t-1)+cx^2(t-1). \end{align*} We obtain $\Phi(\theta)=(1,\theta)$, $\Psi(s)= \begin{bmatrix} -\frac{a_2}{3(\frac{1}{2}a_2-p)^2}-\frac{s}{\frac{1}{2}a_2-p} \\ \frac{1}{\frac{1}{2}a_2-p} \end{bmatrix}$ from Lemma \ref{lem2.4}, and $\frac{1}{2}\tilde{\mathcal{F}}_2(\Phi x,0)=(a+b+c)x_1^2-(b+2c)x_1x_2+cx_2^2$. Applying formulas \eqref{e3.6} and the transforms about \eqref{e3.7}, we deduce that the ODE on center manifold of the origin near $\alpha=0$, \begin{gather*} \dot{z}_1 = z_2, \\ \dot{z}_2 = \rho_1x_1+\rho_2x_2+\eta_1x_1^2+\eta_2x_1x_2. \end{gather*} where \begin{gather*} \rho_1=\frac{1}{\frac{1}{2}a_2-p}\alpha_1 -\frac{2(p-\frac{1}{6}a_2)}{(\frac{1}{2}a_2-p)^2}\alpha_2,\quad \rho_2=-\frac{2(p-\frac{1}{6}a_2)}{(\frac{1}{2}a_2-p)^2}\alpha_1, \\ \eta_1=\frac{a+b+c}{\frac{1}{2}a_2-p}, \quad \eta_2=-\frac{2(p-\frac{1}{6}a_2)(a+b+c)}{(\frac{1}{2}a_2-p)^2} -\frac{b+2c}{\frac{1}{2}a_2-p}. \end{gather*} \end{example} \begin{example}\label{exmp4.2}\rm Consider the following 2-dimensional NFDEs $$\label{e4.2} \begin{gathered} \dot{x}+\dot{x}(t-1) = 2x+y+(1+\alpha_1)x(t-1)+y(t-1)-x^2(t), \\ \dot{y} = \alpha_1x+(1+\alpha_2)y+(\alpha_2-1)y(t-1)-4y^2(t-1). \end{gathered}$$ Comparing with system \eqref{e2.5}, we see that $A(\alpha)= \begin{bmatrix} 2 & 1 \\ \alpha_1 & 1+\alpha_2 \end{bmatrix}, \quad B(\alpha)= \begin{bmatrix} 1+\alpha_1 & 1 \\ 0 & -1+\alpha_2 \end{bmatrix},\quad E= \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}.$ Since $A=A(0)= \begin{bmatrix} 2 & 1\\ 0 &1 \end{bmatrix}$ and $B=B(0)=\begin{bmatrix} 1 & 1\\ 0 &-1 \end{bmatrix}]$, we have $\operatorname{rank}(A+B)=\operatorname{rank} \begin{bmatrix} 3 & 2\\ 0 &0 \end{bmatrix}=1$, condition (i) in Theorem \ref{thm2.1} is valid; on the other hand, by choosing $\phi_1^0= \begin{bmatrix} \frac{1}{3} \\ -\frac{1}{2} \end{bmatrix}$ and $\phi_2^0= \begin{bmatrix} \frac{1}{6} \\ 0 \end{bmatrix}$, one can check that $\operatorname{rank}\Big(A+B, (I+E+B)\phi_1^0\Big) =\operatorname{rank}\begin{bmatrix} 3 & 2 & \frac{1 }{2} \\ 0&0 & 0 \end{bmatrix}=1$ and $\operatorname{rank}\Big(A+B, (I+E+B)\phi_2^0-(E+\frac{1}{2}B)\phi_1^0\Big) =\operatorname{rank} \begin{bmatrix} 3 & 2 & \frac{1 }{4} \\ 0&0 & - \frac{1 }{4} \end{bmatrix}=2$, the conditions (ii'), (iii') in Corollary \ref{cor2.2} are satisfied; then we conclude that system \eqref{e4.2} has a B-T singularity at ($0,0$). Further, we can choose that $\Phi(\theta)= \begin{bmatrix} \frac{1}{3} & \frac{1}{6}+\frac{1}{3}\theta \\ - \frac{1}{2} & - \frac{1}{2}\theta \end{bmatrix}$ and $\Psi(s)= \begin{bmatrix} 0 & -\frac{4}{3}+4s \\ 0 & -4 \end{bmatrix}$ according to Lemma \ref{lem2.4}. By virtue of formulas \eqref{e3.7}, it follows that system \eqref{e4.2} can be reduced to \begin{gather*} \dot{z}_1 = z_2, \\ \dot{z}_2 = (-\frac{4}{3}\alpha_1+4\alpha_2)z_1-(\frac{8}{9}\alpha_1 +\frac{2}{3}\alpha_2)z_2 +\frac{2}{9}z_1^2+\frac{11}{27}z_1z_2. \end{gather*} \end{example} \begin{thebibliography}{00} \bibitem{babram} M. Babram, O. Arino, M. Hbid; Computational scheme of a center manifold for neutral functional differential equations, \emph{J. Math. Anal. Appl.}, 258 (2001) 396-414. \bibitem{chow} S. Chow, C. Li, D. Wang; Normal forms and bifurcation of planar vector fields, Cambridge University Press, \emph{Cambridge}, UK, 2004. \bibitem{faria1} T. Faria, L. T. Magalh\~aes; Normal forms for retarded functional differential equations with parameters and applications to Hopf bifurcation, \emph{J. Differential Equations} 122 (1995) 181-200. \bibitem{faria2} T. Faria, L. T. Magalh\~aes; Normal forms for retarded functional differential equations and applications to Bogdanov-Takens singularity, \emph{J. Differential Equations} 122 (1995) 201-224. \bibitem{guck} J. Guckenheimer, P. Holmes; Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields, \emph{Springer-Verlag}, New York, 1983. \bibitem{hale} J. Hale, S. Lunel; Introduction to functional differential equations, \emph{Springer}, New York, 1993. \bibitem{han} Q. Han; On stability of linear neutral systems with mixed time delays: A discretized Lyapunov functional approach, \emph{Automatica}, 41 (2005) 1209-1218. \bibitem{jiang} W. Jiang, Y. Yuan; Bogdanov-Takens singularity in Van der Pol's oscillator with delayed feedback, \emph{Physica D}, 227 (2007) 149-161. \bibitem{kuzn} Y. Kuznetsov; Elements of applied bifurcation theory, Springer-Verlag, \emph{New York}, 1995. \bibitem{niu1} B. Niu, W. Jiang; Multiple bifurcation analysis in a NDDE arising from van der Pol's equation with extended delay feedback, \emph{Nonlinear Anal. RWA}, 14 (2013) 699-717. \bibitem{niu2} B. Niu, J. Wei; Bifurcation analysis of a NFDE arising from multiple-delay feedback control, \emph{Int. J. Bifurc. Chaos}, 21 (2011) 759-774. \bibitem{qu} Y. Qu, M. Y. Li, J. Wei; Bifurcation analysis in a neutral differential equation, \emph{J. Math. Anal. Appl.}, 378 (2011) 387-402. \bibitem{wang1} C. Wang, J. Wei; Normal forms for NFDEs with parameters and application to the lossless transmission line, \emph{Nonlinear Dyn}, 52 (2008) 199-206. \bibitem{wang2} C. Wang, J. Wei; Hopf bifurcation for neutral functional differential equations, \emph{Nonlinear Anal. RWA}, 11 (2010) 1269-1277. \bibitem{weeder} M. Weedermann; Hopf bifurcation calculati ons for scalar neutral delay differential equations, \emph{Nonlinearity}, 19 (2006) 2091-2102. \bibitem{xu} Y. Xu, M. Huang; Homoclinic orbits and Hopf bifurcations in delay differentianal systems with T-B singularity, \emph{J. Differential Equation}, 244 (2008) 582-598. \bibitem{yuan} Y. Yuan, J. Wei; Singularity analysis on a planar system with multiple delays, \emph{J. Dynam. Diff. Eqs.}, 19(2) (2007) 437-456. \end{thebibliography} \end{document}