\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{amssymb} \AtBeginDocument{{\noindent\small Eighth Mississippi State - UAB Conference on Differential Equations and Computational Simulations. {\em Electronic Journal of Differential Equations}, Conf. 19 (2010), pp. 245--255.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2010 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{245} \title[\hfilneg EJDE-2010/Conf/19/\hfil Bifurcation of solutions] {Bifurcation of solutions of separable parameterized equations into lines} \author[Y.-Q. Shen, T. J. Ypma\hfil EJDE/Conf/19 \hfilneg] {Yun-Qiu Shen, Tjalling J. Ypma} % in alphabetical order \address{Yun-Qiu Shen \newline Department of Mathematics \\ Western Washington University \\ Bellingham, WA 98225-9063, USA} \email{yunqiu.shen@wwu.edu} \address{Tjalling J. Ypma \newline Department of Mathematics \\ Western Washington University \\ Bellingham, WA 98225-9063, USA} \email{tjalling.ypma@wwu.edu} \thanks{Published September 25, 2010.} \subjclass[2000]{65P30, 65H10, 34C23, 37G10} \keywords{Separable parameterized equations; bifurcation; rank deficiency;\hfill\break\indent Golub-Pereyra variable projection method; bordered matrix; singular value decomposition; \hfill\break\indent Newton's method} \begin{abstract} Many applications give rise to separable parameterized equations of the form $A(y, \mu)z+b(y, \mu)=0$, where $y \in \mathbb{R}^n$, $z \in \mathbb{R}^N$ and the parameter $\mu \in \mathbb{R}$; here $A(y, \mu)$ is an $(N+n) \times N$ matrix and $b(y, \mu) \in \mathbb{R}^{N+n}$. Under the assumption that $A(y,\mu)$ has full rank we showed in \cite{shen:ypma:5} that bifurcation points can be located by solving a reduced equation of the form $f(y, \mu)=0$. In this paper we extend that method to the case that $A(y,\mu)$ has rank deficiency one at the bifurcation point. At such a point the solution curve $(y,\mu,z)$ branches into infinitely many additional solutions, which form a straight line. A numerical method for reducing the problem to a smaller space and locating such a bifurcation point is given. Applications to equilibrium solutions of nonlinear ordinary equations and solutions of discretized partial differential equations are provided. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{example}[theorem]{Example} \section{Introduction} Equilibrium solutions of certain nonlinear ordinary differential equations with parameters, and numerical solutions of certain boundary value problems, give rise to separable nonlinear parameterized equations of the form $$A(y, \mu)z+b(y, \mu)=0 \label{eq:1.1}$$ where $y \in {\mathbb{R}}^n, z \in {\mathbb{R}}^N$ and the parameter $\mu \in {\mathbb{R}}$. Here $A(y, \mu)$ is an $(N+n) \times N$ matrix and $b(y, \mu) \in {\mathbb{R}}^{N+n}$; in practice usually $N \gg n$. We refer to the recent survey \cite{golu:pere:2} for numerous examples, references, and variants of the Golub-Pereyra variable projection method for the solution of separable equations. Further references to related material appear in \cite{golu:leve, golu:pere:1,luke, morr, mull:stok,shen:ypma,ypma:shen}. Recently in \cite{shen:ypma:5} we extended a variant of the Golub-Pereyra variable projection method \cite{shen:ypma,ypma:shen} to \eqref{eq:1.1} for bifurcation problems specified by a parameter $\mu$. Under the assumption that $A(y,\mu)$ has full rank we show there that curve following and bifurcation at a point $(y, \mu,z)$ can be located within a smaller space, for $(y, \mu)$ only, by solving a reduced equation of the form $f(y, \mu) = 0$. In this paper we extend both the technique and the analysis of \cite{shen:ypma:5} to the case when $A(y, \mu)$ has rank deficiency one at the bifurcation point $(y^*,\mu^*,z^*)$. That implies a line of solutions in the nullspace of $A(y^*,\mu^*)$ satisfying \eqref{eq:1.1} through the bifurcation point. We call this phenomenon bifurcation into a line. We show that \eqref{eq:1.1} can be reduced to solving only for $(y, \mu, w)$, for a variable $w \in \mathbb{R}$ defined below, by using a variant of the matrix bordering techniques of \cite{beyn,grie:redd,rabi:redd,shen,shen:2,shen:3,shen:ypma:3} and an extension of the algorithms of \cite{shen:ypma:2,shen:ypma:4} to eliminate $N-1$ components of the linear variable $z$. We show that bifurcation phenomena, including the line of solutions, are preserved in this smaller space and that the relevant points can be located by solving a smaller separable nonlinear system of the form $$f(y,\mu,w) \equiv \alpha(y,\mu)w+ \beta(y,\mu) = 0, \label{eq:1.2}$$ where $w \in \mathbb{R}$ and $\alpha(y,\mu), \beta(y,\mu) \in {\mathbb{R}}^{n+1}$. If $N$ is much larger than $n$ the reduction of dimension from $N+n+1$ to $n+2$ is significant. The precise form of $f(y,\mu,w)$ and how to solve \eqref{eq:1.2} are detailed below. In the next section we present the theoretical basis of our method. A reduction from \eqref{eq:1.1} to \eqref{eq:1.2} is presented and the preservation of the bifurcation is proved. In Section 3 we derive expressions for the derivatives involved in our method. We derive an extended system for the location of the bifurcation point in the smaller dimensional space in Section 4. In the last section we give some numerical examples to illustrate the method. \section{Analysis} Write $$\widetilde{f} (y, \mu, z) \equiv A(y, \mu)z+b(y, \mu)=0 \label{eq:2.1}$$ where $y \in {\mathbb{R}}^n, z \in {\mathbb{R}}^N, \mu \in {\mathbb{R}}$ and $A(y, \mu)$ is an $(N+n) \times N$ matrix and $b(y, \mu) \in {\mathbb{R}}^{N+n}$. Assume both $A$ and $b$ are $C^2-Lipschitzian$ in $y$ and $\mu$; thus $\widetilde{f}(y,\mu,z)$ is also $C^2$-Lipschitzian in $y$ and $\mu$. Suppose a curve following technique leads to a bifurcation point $(y^*, \mu^*,z^*)$ of the solution set where the column rank deficiency of $A(y^*, \mu^*)$ is one. At this point bifurcation into a straight line in the nullspace of $A(y^*,\mu^*)$ and other curves occurs. The singular value decomposition of $A(y, \mu)$ at any point $(y, \mu)$ is defined to be \cite{golu:vanl,watk} $$A(y, \mu)=U \Sigma V^T \equiv [u_1,\dots,u_{N+n}] \begin{bmatrix} \operatorname{diag}(\sigma_1, \dots, \sigma_N) \\ 0_{n \times N} \end{bmatrix} [v_1,\dots,v_N]^T \label{eq:2.2}$$ where $\{ u_1,\dots,u_{N+n} \}$ and $\{ v_1,\dots,v_N \}$ form orthonormal bases for ${\mathbb{R}}^{N+n}$ and ${\mathbb{R}}^N$ respectively and the singular values satisfy $\sigma_1 \ge \dots \ge \sigma_N \ge 0$. The rank deficiency one at $(y^*,\mu^*)$ implies that the corresponding singular values satisfy $\sigma^*_{N-1} \ne 0$ and $\sigma^*_N=0$. If $(\bar{y}, \bar{\mu})$ is near $(y^*,\mu^*)$ then the corresponding singular values satisfy $\bar{\sigma}_{N-1} \ne 0$ and $\bar{\sigma}_N$ is small by continuity of the singular values. For a given (fixed) $(\bar y, \bar{\mu})$ near $(y^*, \mu^*)$ define the $N$-vector $\mathit{l}$ and the $(N+n) \times (n+1)$ matrix $R$ from the singular value decomposition of $A(\bar{y}, \bar{\mu})=\bar{U} \bar{\Sigma} {\bar V}^T$ by $$\mathit{l} \equiv \bar{v}_N, \quad R \equiv [\bar{u}_N, \dots, \bar{u}_{N+n}]. \label{eq:2.3}$$ Let $M(y,\mu)$ be the $(N+n+1) \times (N+n+1)$ bordered matrix $$M(y, \mu) \equiv \begin{bmatrix} A(y, \mu) & R \\ \mathit{l}^T & O \end{bmatrix}. \label{eq:2.4}$$ Note that in \eqref{eq:2.4} and henceforth the vector $\mathit{l}$ and the matrix $R$ are fixed; they are determined by the selected point $(\bar y, \bar\mu)$; only the term $A(y, \mu)$ varies with $y$ and $\mu$. \begin{lemma}\label{lm:2.1} The matrix $M(y, \mu)$ is nonsingular for all $(y, \mu)$ in a small neighborhood of $(y^*, \mu^*)$. \end{lemma} \begin{proof} Let $\bar{y}, \bar{\mu}$ be as above. Observe that since $\bar{\sigma}_{N-1} \ne 0$ $M(\bar y, \bar{\mu}) = \begin{bmatrix} \bar{U} & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \operatorname{diag}(\bar{\sigma}_1, \dots, \bar{\sigma}_{N-1}) & 0 & 0 \\ 0 & \begin{bmatrix} \bar{\sigma}_N \\ 0 \end{bmatrix} & I_{n+1} \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \bar{V}^T & 0 \\ 0 & I_{n+1} \end{bmatrix}$ is nonsingular. From the Neumann Lemma \cite{orte:rhei} $M(y, \mu)$ is nonsingular provided $\| M(\bar y, \bar {\mu})^{-1}[M(\bar y, \bar {\mu})-M(y, \mu)]\| < 1$, and since $\| M(\bar y, \bar {\mu})-M(y, \mu) \| = \Big\| \begin{bmatrix} A(\bar y, \bar {\mu})-A(y, \mu) & 0 \\ 0 & 0 \end{bmatrix} \Big\|$ this holds for all $(y, \mu)$ and $(\bar y, \bar{\mu})$ sufficiently near $(y^*, \mu^*)$ by continuity of $A(y, \mu)$. \end{proof} Assume $z$ satisfies \eqref{eq:1.1}. Define $w=\mathit{l}^Tz \in \mathbb{R}$ and $\tilde z=z-w \mathit{l}$. Then $$\begin{bmatrix} A(y,\mu) \\ \mathit{l}^T \end{bmatrix} z + \begin{bmatrix} b(y,\mu) \\ -w \end{bmatrix} =\begin{bmatrix}0 \\ 0 \end{bmatrix} =\begin{bmatrix} A(y,\mu) \\ \mathit{l}^T \end{bmatrix} (\tilde z + w\mathit{l})+ \begin{bmatrix} b(y,\mu) \\ -w \end{bmatrix}. \label{eq:2.5}$$ Denote $$\widetilde A(y, \mu)=\begin{bmatrix} A(y, \mu) \\ \mathit{l}^T \end{bmatrix}, \quad \widetilde b(y, \mu,w) =\begin{bmatrix} b(y, \mu) \\ -w \end{bmatrix}. \label{eq:2.6}$$ By Lemma \ref{lm:2.1} the first $N$ columns of $M(y,\mu)$ form the full rank matrix $\widetilde A(y, \mu)$. Thus \eqref{eq:2.5} has the unique least squares solution $\widetilde{z}+w\mathit{l}=-[\widetilde A(y, \mu)]^+\widetilde{b}(y, \mu,w)$, where $[\widetilde{A}(y, \mu)]^+$ denotes the pseudo-inverse of $\widetilde A(y, \mu)$. Extending the algorithm of \cite{shen:ypma,ypma:shen} we have the following result. \begin{theorem}\label{th:2.1} $(y, \mu, w, \tilde z)$ is a solution of \eqref{eq:2.5} if and only if it satisfies $$f(y,\mu,w) \equiv C^T(y, \mu) \widetilde{b}(y,\mu,w)= 0, \quad \tilde z = -w\mathit{l}- [\widetilde{A}(y, \mu)]^+ \widetilde{b}(y,\mu,w), \label{eq:2.7}$$ where $C^T(y,\mu)$ is the transpose of the $(N+n+1) \times (n+1)$ matrix $C(y,\mu)$ satisfying $$C^T(y, \mu)M(y, \mu) = C^T(y, \mu) \begin{bmatrix} A(y, \mu) & R \\ \mathit{l}^T & 0 \end{bmatrix} = -[0_{(n+1) \times N}, I_{n+1}]. \label{eq:2.8}$$ \end{theorem} \begin{proof} ($\Longrightarrow$) Assume \eqref{eq:2.5} holds. Multiply both sides of \eqref{eq:2.5} by $C^T(y, \mu)$ and use \eqref{eq:2.8} to give the first part of \eqref{eq:2.7}. Since $z=\tilde z +w\mathit{l}$ is a solution of \eqref{eq:2.5} for the given $(y,\mu,w)$ it is the unique least squares solution of a linear equation with a full rank coefficient matrix, so we have the second part of \eqref{eq:2.7}. ($\Longleftarrow$) Assume \eqref{eq:2.7} holds. From \eqref{eq:2.8} we have $0= C^T(y, \mu) \widetilde b(y, \mu,w) =-[0, I_{n+1}]M^{-1}(y, \mu) \widetilde{b}(y, \mu,w)$ which implies \begin{aligned} \widetilde{b}(y, \mu,w) &=M(y, \mu) \begin{bmatrix} I_N & 0 \\ 0 & I_{n+1} \end{bmatrix} M^{-1}(y,\mu) \widetilde{b}(y, \mu,w) \\ &=M(y, \mu) \begin{bmatrix} [I_N, 0]M^{-1}(y,\mu) \widetilde{b}(y, \mu,w) \\ 0 \end{bmatrix} \\ &= \widetilde{A}(y, \mu) [I_N, 0]M^{-1}(y,\mu) \widetilde{b} (y, \mu,w). \end{aligned} \label{eq:2.9} We then obtain, using \eqref{eq:2.7} and \eqref{eq:2.9}, \begin{align*} &\widetilde A(y, \mu)(\tilde{z}+w\mathit{l}) + \widetilde b(y, \mu,w) \\ &= \widetilde A(y, \mu) \left( -{\widetilde A}^+ (y, \mu)\widetilde b(y, \mu,w) \right) + \widetilde{b}(y, \mu,w) \\ &= [I_{N+n+1}-\widetilde{A}(y, \mu)\widetilde{A}^+(y, \mu)]\widetilde{A}(y, \mu) [I_N, 0]M^{-1}(y,\mu)\widetilde{b}(y, \mu,w) \\ &= [\widetilde{A}(y, \mu)- \widetilde{A}(y, \mu)\widetilde{A}^+ (y, \mu)\widetilde{A}(y, \mu)][I_N, 0]M^{-1}(y,\mu)\widetilde{b} (y, \mu,w) = 0 \end{align*} as claimed. \end{proof} We only need to solve the first equation of \eqref{eq:2.7} for $(y,\mu,w)$ since it preserves bifurcations of \eqref{eq:2.5} by the following theorem. \begin{theorem} \label{th:2.2} For any fixed parameter $\mu$, the points $(y, \mu,w)$ and $(\bar{y},\mu, \bar{w})$ are two distinct solutions of the first part of \eqref{eq:2.7} if and only if $(y, \mu,z)$ and $(\bar y, \mu, \bar z)$ are two distinct solutions of \eqref{eq:2.5} with $z=\tilde{z}+w \mathit{l}$ and $\bar z = \tilde{\bar{z}}+ \bar{w}\mathit{l}$, where $\tilde{z}$ and $\tilde{\bar{z}}$ are determined by the second part of \eqref{eq:2.7} respectively. \end{theorem} \begin{proof} $( \Longrightarrow )$ If not, then $(y,\mu,w) \ne (\bar{y},\mu, \bar{w})$ but $(y, \mu,z)=( \bar{y}, \mu, \bar{z})$ which implies that $y=\bar{y}$, $z= \bar{z}$ and thus $w \ne \bar{w}$. Thus $0=z-\bar{z}=(\tilde{z}-\tilde{\bar z})+(w-\bar{w})\mathit{l}$ by the definition of $z$ and $\bar{z}$. From Theorem \ref{th:2.1}, the fact that $(y,\mu,z)$ satisfies the first equation of \eqref{eq:2.7} and $\tilde z$ satisfies the second equation of \eqref{eq:2.7} implies that $(y,\mu,w,\tilde{z})$ satisfies \eqref{eq:2.5}; similarly so does $(\bar{y},\mu,\bar{w}, \tilde{\bar z})$. Therefore, $$\widetilde{A}(y, \mu)(\tilde z + w\mathit{l}) + \widetilde{b}(y,\mu,w)=0 =\widetilde{A}(\bar{y}, \mu)(\tilde{\bar z} + \bar{w}\mathit{l}) +\widetilde{b}(\bar{y}, \mu, \bar{w}).$$ Subtracting the left side of the equation from the right side and using $y = \bar{y}$ and $z = \bar{z}$ we have $w=\bar{w}$, contradicting $w \ne \bar{w}$. So $(y,\mu,z) \ne ( \bar{y}, \mu, \bar{z})$. $( \Longleftarrow )$ If not, then $(y,\mu,z) \ne ( \bar{y}, \mu, \bar{z})$ but $(y,\mu,w) = (\bar{y},\mu, \bar{w})$ which implies that $w= \bar{w}$, $y= \bar{y}$ but $z \ne \bar{z}$. By the second equation of \eqref{eq:2.7} $\tilde{z} = \tilde{\bar{z}}$. Thus $z=\tilde{z}+w\mathit{l}=\tilde{\bar{z}} +\bar{w}\mathit{l}=\bar{z}$, which is a contradiction. So $(y,\mu,w) \ne (\bar{y},\mu, \bar{w})$. \end{proof} Write $$C^T(y,\mu) \equiv [(P^T(y, \mu))_{(n+1) \times (N+n)},(q^T(y,\mu))_{n+1}]. \label{eq:2.10}$$ Then the reduced equation $f(y,\mu,w)=0$ of \eqref{eq:2.7} becomes $$f(y,\mu,w) \equiv - q^T(y,\mu)w+ P^T(y,\mu) b(y,\mu) \label{eq:2.11}$$ which is of the form of a separable equation in a smaller space. Essentially, the $N$-vector $z$ has been has been replaced by the scalar $w$ while the matrix $(A(y,\mu))_{(N+n) \times N}$ is replaced by the $(n+1)$-vector $-q^T(y,\mu)$. Finally we have the following results showing preservation of the line of solutions. \begin{theorem} \label{th:2.3} Let $(y^*, \mu^*, z^*)$ be a bifurcation point of \eqref{eq:2.1} with $A(y^*,\mu^*)$ having rank deficiency one and let $A(y,\mu)$ have full rank at every solution point $(y,\mu,z)$ of \eqref{eq:2.1} in a small neighborhood of $(y^*, \mu^*, z^*)$ with $(y,\mu) \ne (y^*, \mu^*)$. Let $(\hat{y},\hat{\mu},\hat{z})$ be a solution point of \eqref{eq:2.1} in this neighborhood. Then $(\hat{y},\hat{\mu})=(y^*,\mu^*)$ if and only if $q^T(\hat{y},\hat{\mu})=0 \in \mathbb{R}^{n+1}$. \end{theorem} \begin{proof} First note from \eqref{eq:2.8} that $C(y,\mu)$ has full rank $n+1$ since $M(y,\mu)$ is nonsingular for all $(y,\mu)$ in a neighborhood of $(y^*,\mu^*)$. $( \Longrightarrow )$ From \eqref{eq:2.8} the rows of $C^T(y^*,\mu^*)=[P^T(y^*,\mu^*),q^T(y^*,\mu^*)]$ form a basis for the left nullspace of $\widetilde A(y^*, \mu^*)$, which is $(n+1)$-dimensional because the matrix $M(y^*,\mu^*)$ of \eqref{eq:2.4} has full rank. Since $A(y^*,\mu^*)$ has column rank deficiency one, the matrix $A(y^*,\mu^*)$ also has an $(n+1)$-dimensional left nullspace. Suppose its basis consists of the vectors $\varphi_1,\dots, \varphi_{n+1}$, and denote by $\Phi$ the matrix with these vectors as columns. Then the rows of $[ \Phi^T,0_{(n+1) \times 1}]$ also form a basis for the left nullspace of $\widetilde{A}(y^*, \mu^*)$, so $[P^T(y^*,\mu^*),q^T(y^*,\mu^*)]=G [ \Phi^T,0]$ for some nonsingular $(n+1) \times (n+1)$ matrix $G$, which implies that $$P^T(y^*,\mu^*)=G \Phi^T, \quad q^T(y^*,\mu^*)=0.$$ $( \Longleftarrow )$ The $n+1$ rows of $C^T(\hat{y},\hat{\mu})=[P^T(\hat{y},\hat{\mu}),0]$ form a basis of the left nullspace of $\widetilde A(\hat{y},\hat{\mu})$, which implies that the $n+1$ rows of $P^T(\hat{y},\hat{\mu})$ form a basis for the left nullspace of the $(N+n) \times N$ matrix $A(\hat{y},\hat{\mu})$. Thus $A(\hat{y},\hat{\mu})$ has rank deficiency one which implies that $(\hat{y},\hat{\mu})=(y^*,\mu^*)$. \end{proof} \begin{theorem} \label{th:2.4} (a) If $(y^*,\mu^*,z^*)$ is a bifurcation point of \eqref{eq:1.1} which bifurcates solutions into a straight line then $(y^*,\mu^*,w^*)$ with $w^*=\mathit{l}^Tz^*$ is a bifurcation point of \eqref{eq:2.11} which bifurcates solutions into a straight line; (b) If $(y^*,\mu^*,w^*)$ is a bifurcation point of \eqref{eq:2.11} which bifurcates solutions into a straight line then $(y^*,\mu^*,z^*)$ with $z^*=-[\widetilde{A}(y^*, \mu^*)]^+ \widetilde{b}(y^*,\mu^*,w^*)$ is a bifurcation point of \eqref{eq:1.1} which bifurcates solutions into a straight line. \end{theorem} \begin{proof} From Theorems \ref{th:2.1} and \ref{th:2.2} $(y^*,\mu^*,z^*)$ is a bifurcation point of \eqref{eq:1.1} if only if $(y^*,\mu^*,w^*)$ is a bifurcation point of \eqref{eq:2.11}. From Theorem \ref{th:2.3} there is a line of solutions of \eqref{eq:1.1} through $(y^*,\mu^*,z^*)$ if and only if there is a line of solutions of \eqref{eq:2.11} through $(y^*,\mu^*,w^*)$. \end{proof} \section{Computation} \label{sc:3} We shall need derivatives of $f(y,\mu,w)$ to solve \eqref{eq:2.11}. Denote the vectors $f'_j \equiv \frac{\partial f}{\partial {y_j}}$ for $j=1,\dots,n$, $f'_{n+1} \equiv \frac{\partial f}{\partial \mu}$ and $f'_{n+2} \equiv \frac{\partial f}{\partial w}$. We use similar notation for the vectors of second partial derivatives $f''_{j,k}=(f'_j)'_k$ of $f$ with respect to the $j$th and $k$th independent variables of $(y,\mu,w)$, and also for the derivatives of other functions of $(y,\mu,w)$. \begin{theorem} \label{th:3.1} Let $\zeta(y,\mu,w) \in \mathbb{R}^N$ and $\psi(y,\mu,w) \in \mathbb{R}^{n+1}$ satisfy $$\begin{bmatrix} A(y, \mu) & R \\ \mathit{l}^T & 0 \end{bmatrix} \begin{bmatrix} \zeta(y,\mu,w) \\ \psi(y,\mu,w) \end{bmatrix} = - \begin{bmatrix} b(y,\mu) \\ -w \end{bmatrix}. \label{eq:3.1}$$ Then for $j,k=1,\dots,n+1$, we have \begin{itemize} \item[(a)] $$f(y,\mu,w)= \psi(y,\mu,w) \label{eq:3.2}$$ \item[(b)] $f'_j=\psi'_j = P^T(y, \mu) [ A'_j \zeta(y,\mu,w)+b'_j]$, $$\zeta'_j=-[I_N,0_{N \times (n+1)}]M^{-1}(y,\mu) \begin{bmatrix} A'_j \zeta(y,\mu,w)+b'_j \\ 0 \end{bmatrix}; \label{eq:3.3}$$ \item[(c)] $f'_{n+2}=\psi'_{n+2}=-q^T(y,\mu)$, $$\zeta'_{n+2} = [I_N,0] M^{-1}(y, \mu) e_{N+n+1}; \label{eq:3.4}$$ \item[(d)] $$f''_{j,k}= \psi''_{j,k} = P^T(y,\mu)[A''_{j,k}\zeta(y,\mu,w) +A^{ \prime}_j \zeta'_{k} +A'_k \zeta'_{j}+b''_{j,k}]; \label{eq:3.5}$$ \item[(e)] $f''_{n+2,j}= \psi''_{n+2,j} = P^T(y, \mu) A'_j \zeta'_{n+2}$, $$f''_{j,n+2}=\psi''_{j,n+2} = P^T(y, \mu) A'_j \zeta'_{n+2}; \label{eq:3.6}$$ \item[(f)] $$f''_{n+2,n+2}=\psi''_{n+2,n+2} = 0 \in \mathbb{R}^{n+1}. \label{eq:3.7}$$ \end{itemize} where $e_{N+n+1}$ is the $(N+n+1)$-th unit vector of the standard basis in $\mathbb {R}^{N+n+1}$. \end{theorem} \begin{proof} (a) Using \eqref{eq:2.7}, \eqref{eq:2.8} and \eqref{eq:3.1} we obtain $$f(y,\mu,w)=C^T(y,\mu)\widetilde{b}(y,\mu,w)= -[0,I_{n+1}]M^{-1}(y,\mu)\widetilde{b}(y,\mu,w) =\psi(y,\mu,w).$$ (b) Differentiating both sides of \eqref{eq:3.1} with respect to $y_j$ or $\mu$, we have $$\begin{bmatrix} A'_j & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} \zeta(y,\mu,w) \\ \psi(y,\mu,w) \end{bmatrix} + M(y,\mu) \begin{bmatrix} \zeta'_j \\ \psi'_j \end{bmatrix} = - \begin{bmatrix} b'_j \\ 0 \end{bmatrix} \label{eq:3.8}$$ which implies $\begin{bmatrix} \zeta'_j \\ \psi'_j \end{bmatrix} = -M^{-1}(y,\mu) \begin{bmatrix} A'_j \zeta(y,\mu,w)+b'_j \\ 0 \end{bmatrix}.$ Using \eqref{eq:2.8}, \eqref{eq:2.10} and \eqref{eq:3.2} we thus have \begin{align*} f'_j &= \psi'_j =-[0, I_{n+1}]M^{-1}(y,\mu) \begin{bmatrix} A'_j \zeta(y,\mu,w)+b'_j \\ 0 \end{bmatrix} \\ &= P^T(y, \mu) [A'_j \zeta(y,\mu,w)+b'_j] \end{align*} and $\zeta'_j=-[I_N,0]M^{-1}(y,\mu) \begin{bmatrix} A'_j \zeta(y,\mu,w)+b'_j \\ 0 \end{bmatrix}.$ (c) Differentiating both sides of \eqref{eq:3.1} with respect to $w$ gives $$M(y, \mu) \begin{bmatrix} \zeta'_{n+2} \\ \psi'_{n+2} \end{bmatrix} = e_{N+n+1} \label{eq:3.9}$$ which with \eqref{eq:2.8}, \eqref{eq:2.10} and \eqref{eq:3.2} gives $f'_{n+2}=\psi'_{n+2} = [0,I_{n+1}] M^{-1}(y, \mu) e_{N+n+1} = -C^T(y, \mu)e_{N+n+1} = -q^T(y,\mu)$ and $\zeta'_{n+2} = [I_N,0] M^{-1}(y, \mu) e_{N+n+1}$. (d) Differentiating both sides of \eqref{eq:3.8} with respect to $y_k$ or $\mu$ with \eqref{eq:2.8}, \eqref{eq:2.10} and \eqref{eq:3.2} we obtain \eqref{eq:3.5}. (e) Differentiating both sides of \eqref{eq:3.9} with respect to $y_j$ or $\mu$ with \eqref{eq:2.8}, \eqref{eq:2.10} and \eqref{eq:3.2} we obtain the first part of \eqref{eq:3.6}. Direct differentiation of both sides of the first part of \eqref{eq:3.3} with respect to $w$ gives the second part of \eqref{eq:3.6} (f) Differentiating both sides of the first part of \eqref{eq:3.4} with respect to $w$ we obtain \eqref{eq:3.7}. \end{proof} Let $(y^*,\mu^*,w^*)$ be a bifurcation solution of \eqref{eq:2.11}. By \eqref{eq:3.2} $\psi(y^*,\mu^*,w^*) =f(y^*,\mu^*,w^*)=0$. Thus \eqref{eq:2.5} and \eqref{eq:3.1} imply that \begin{aligned} z^* &= -[\widetilde{A}(y^*, \mu^*)]^+ \widetilde{b}(y^*,\mu^*,w^*) \\ &= \zeta(y^*,\mu^*,w^*)\\ & = -[I_N,0_{N \times (n+1)}]M^{-1}(y^*,\mu^*) \widetilde{b}(y^*,\mu^*,w^*) \qquad \end{aligned}\label{eq:3.10} by the uniqueness of the linear least squares solution for the full rank matrix $\widetilde{A}(y^*,\mu^*)$. This expression is used to obtain $z^*$, without using the pseudo-inverse, after $y^*,\mu^*$ and $w^*$ have been found. From Theorem \ref{th:3.1} we obtain the $(n+1) \times (n+2)$ Jacobian matrix $$f'(y,\mu,w) =[P^T(y,\mu)K(y,\mu,w),-q^T(y,\mu)] = C^T(y, \mu) \begin{bmatrix} K(y,\mu,w)& 0 \\ 0 & -1 \end{bmatrix} \label{eq:3.11}$$ where the $(N+n) \times (n+1)$ matrix $$K(y,\mu,w) \equiv [A'_1 \zeta(y,\mu,w), \dots, A'_{n+1}\zeta(y,\mu,w)]+b'(y,\mu) \label{eq:3.12}$$ and $P(y, \mu)$ and $q(y,\mu)$ are defined in \eqref{eq:2.11}. Since $(y^*,\mu^*,w^*)$ is a bifurcation point of \eqref{eq:2.11}, $f'(y^*,\mu^*,w^*)$ has a rank deficiency $d \ge 1$. To solve the equation \eqref{eq:2.11} we use the technique of \cite{shen:2,shen:3,shen:ypma:3,shen:ypma:5}. Let the singular value decomposition of the $(n+1) \times (n+2)$ matrix $f'(\bar{y},\bar{\mu},\bar{w})$ be \begin{aligned} & f'(\bar{y},\bar{\mu},\bar{w}) =\widetilde{U} \widetilde{\Sigma} {\widetilde V}^T \\ & \equiv [\widetilde{u}_1,\dots,\widetilde{u}_{n+1}] [ \operatorname{diag}(\widetilde{\sigma}_1, \dots, \widetilde {\sigma}_{n+1}) , 0_{(n+1) \times 1}] [\widetilde{v}_1,\dots,\widetilde{v}_{n+2}]^T. \end{aligned}\label{eq:3.13} Let the $(n+2) \times (d+1)$ matrix $\widetilde{L}$ and the $(n+1) \times d$ matrix $\widetilde{R}$ be $$\widetilde{L} = [\widetilde{v}_{n+2-d},\dots,\widetilde{v}_{n+2}], \quad \widetilde{R} =[\widetilde {u}_{n+2-d},\dots,\widetilde{u}_{n+1}]. \label{eq:3.14}$$ An argument similar to that used in Lemma \ref{lm:2.1} leads to \begin{lemma}\label{lm:3.1} The $(n+d+2) \times (n+d+2)$ bordered matrix $$\widetilde{M}(y, \mu,w) \equiv \begin{bmatrix} f'(y,\mu,w) & \widetilde{R} \\ \widetilde{L}^T & 0 \end{bmatrix} \label{eq:3.15}$$ is nonsingular for all $(y, \mu,w)$ in a small neighborhood of $(y^*, \mu^*,w^*)$. \end{lemma} As in \cite{shen:3,shen:ypma:3} an extended system of $n+d+2$ equations in $n+d+2$ variables for \eqref{eq:2.11} is $$F(y,\mu,w,\lambda)= \begin{bmatrix} f(y,\mu,w)+ \widetilde{R} \lambda \\ g(y,\mu,w) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} , \label{eq:3.16}$$ with the Jacobian matrix $$F'(y,\mu,w,\lambda)= \begin{bmatrix} f'(y,\mu,w) & \widetilde{R}\\ -\eta^T(y,\mu,w) \sum_{k=1}^{n+1} [\xi_k(y,\mu,w) \nabla^2f_k(y,\mu,w)] & 0 \end{bmatrix}. \label{eq:3.17}$$ Here $\lambda \in \mathbb{R}^{d}$ is a Lagrange-type auxiliary vector, $\xi(y,\mu,w) \in \mathbb{R}^{n+1}$ and $g(y,\mu,w) \in \mathbb{R}^{d+1}$ satisfy $$\widetilde{M}^T(y,\mu,w) \begin{bmatrix} \xi(y,\mu,w) \\ g(y,\mu,w) \end{bmatrix} = \begin{bmatrix} 0 \\ \gamma \end{bmatrix} \label{eq:3.18}$$ with a randomly selected nonzero vector $\gamma \in \mathbb{R}^d$, and $\eta(y,\mu,w) \in \mathbb{R}^{(n+1) \times (d+1)}$ satisfies $$\eta(y,\mu,w)= \begin{bmatrix} I_{n+1} & 0 \end{bmatrix} \widetilde{M}^{-1}(y,\mu,w) \begin{bmatrix} 0 \\ I_{d+1} \end{bmatrix}. \label{eq:3.19}$$ It has been proved that under certain conditions for almost every choice of $\gamma$ the matrix $F'(y^*,\mu^*,w^*,0)$ is nonsingular (see \cite{shen:3,shen:ypma:3} for details). Thus Newton's method can be used to solve \eqref{eq:3.16} with quadratic convergence. The solution of \eqref{eq:2.11} corresponds to a solution of \eqref{eq:3.16} with $\lambda^*=0$. Convergence of Newton's method depends on the initial point being near the solution. We can choose a good initial point by using a curve following technique on a solution branch of \eqref{eq:2.11} connecting to the bifurcation point \cite{chow:shen,shen:ypma:5}. That starting point corresponds to a starting point of \eqref{eq:3.16} with $\lambda^{(0)}=0$. The correct rank deficiency $d$ can be obtained by looking for small singular values $\widetilde{\sigma}_{n+2-d},\dots,\widetilde{\sigma}_{n+1}$ in \eqref{eq:3.13}, since the singular values are continuous and the corresponding values at $(y^*,\mu^*,w^*)$ are zero. If the rank deficiency $d$ is under-estimated then convergence becomes linear, while if the rank deficiency $d$ is over-estimated then the computed solution may not be the desired solution with $\lambda^*=0$. \section{Examples} We present two examples to illustrate our method. \begin{example}\label{ex:1} \rm $$\frac{d}{dt} \begin{bmatrix} z_1 \\ y \\ z_2 \end{bmatrix} = \begin{bmatrix} -\mu & 0 \\ -1 & \mu y \\ 0 & - \mu \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} + \begin{bmatrix} y \\ -\mu y \\ \mu -\mu y^2 \end{bmatrix}. \label{eq:5.1}$$ \end{example} This simple example with $N=2$ and $n=1$ is a case of the Sherman system \cite{sher} with one parameter $\mu$ while the other two parameters are fixed. The system is derived from a nuclear spin generator; for more details see \cite{niko:bozh:nede:zlat}. There is one equilibrium solution $y=0,z_1=0, z_2=1$ for any value of $\mu$. At the point $(\bar y,\bar \mu,\bar z_1,\bar z_2)=(0.02,0.02,0.02,1.02,)$ near the bifurcation point $(y^*, \mu^*,z_1^*, z_2^*)=(0,0,0,1)$ the singular value decomposition of $A(\bar y, \bar{\mu})$ is $$\begin{bmatrix} -0.0200 & 0.0004 &-0.9998 \\ -0.9998 & 0.0000 & 0.0200 \\ 0.0000 & 1.0000 & 0.0004 \end{bmatrix} \begin{bmatrix} 1.0002 & 0 \\ 0 & 0.0200 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1.0000 & -0.0004 \\ -0.0004 & -1.0000 \end{bmatrix}^T.$$ $A(y^*,\mu^*)$ has rank deficiency one as indicated by the singular values of $A(\bar{y}, \bar{\mu})$ since $\bar{\sigma}_2=0.0200$ is very small relative to $\bar{\sigma}_1=1.0002$. This gives $$\mathit{l}=\begin{bmatrix} -0.0004 \\ -1.0000 \end{bmatrix}, \quad R= \begin{bmatrix} 0.0004 & -0.9998 \\ 0.0000 & 0.0200 \\ 1.0000 & 0.0004 \end{bmatrix}$$ and $\bar w=\mathit{l}^T \bar z=-1.0200$. The singular value decomposition of $f'(\bar y,\bar \mu, \bar w)$ is $\begin{bmatrix} 0.0004 & 1.0000 \\ 1.0000 & -0.0004 \end{bmatrix} \begin{bmatrix} 0.9998 & 0 & 0 \\ 0 & 0.0286 & 0 \end{bmatrix} \begin{bmatrix} 1.0000 & 0.0003 & 0.0003 \\ -0.0004 & 0.7142 & 0.6999 \\ 0.0000 & -0.6999 & 0.7142 \end{bmatrix}^T$ which suggests that the rank deficiency of $f'(y^*,\mu^*,w^*)$ is $d=1$ since $\widetilde \sigma_2=0.0286$ is very small relative to $\widetilde \sigma_1=0.9998$. Hence $\widetilde L=\begin{bmatrix} 0.0003 & 0.003 \\ 0.7142 & 0.6999 \\ -0.6999 & 0.7142 \end{bmatrix}, \quad \widetilde R = \begin{bmatrix} 1.0000 \\ -0.0004 \end{bmatrix}.$ We choose $\gamma=0.6742$ at random. We list in Table \ref{tb:1} the data beginning with $(y^{(0)},\mu^{(0)},w^{(0)},\lambda^{(0)})=(0.02,0.02,-1.02,0)$ and using $$\| d^{(k)}\|_2 = \|(y^{(k+1)},\mu^{(k+1)},w^{(k+1)},\lambda^{(k+1)})-(y^{(k)}, \mu^{(k)},w^{(k)}, \lambda^{(k)})\|_2.$$ \begin{table}[thbp] \caption{Newton iterations for computing a bifurcation point in Example \ref{ex:1}} \label{tb:1} \begin{center} \begin{tabular}{||c||c|c|c|c||c||} \hline\hline $k$ & $y^{(k)}$ & $\mu^{(k)}$ & $w^{(k)}$ & $\lambda^{(k)}$ & $\|d^{(k)}\|_2$ \\ \hline\hline 0 & ~2.0000e-02 & ~2.0000e-02 & -1.0200e+00 & ~0.0000e+00 & 3.4423e-02 \\ 1 & -8.1505e-06 & -6.2170e-09 & -1.0004e+00 & ~4.0814e-04 & 5.7221e-04 \\ 2 & -2.9878e-15 & -8.4432e-19 & -1.0000e+00 & -2.4929e-12 & 1.3296e-10 \\ 3 & ~0.0000e+00 & -4.7698e-26 & -1.0000e+00 & ~0.0000e+00 & ------------ \\ \hline\hline \end{tabular} \end{center} \end{table} The last line of the table is an approximation of $(y^*,\mu^*,w^*,\lambda^*)$; using \eqref{eq:3.10} then gives $z^* \approx (0.0000e+00,1.0000e+00)^T$. \begin{example} \label{ex:2} \rm Let $(u,v)=(u(t,x),v(t))$ be a solution of the system $$\begin{gathered} \frac{\partial u}{\partial t}=\mu u + \frac{\partial^2u}{\partial x^2} -c_1uv, \quad \frac{dv}{dt}=(\mu-\mu_0)v-c_2v^2+c_3v \int_{0}^1 u(t,x)dx, \\ x \in (0,1), \quad u(t,0)=u(t,1)=0 \label{eq:4.2} \end{gathered}$$ where $c_1,c_2$ and $c_3$ are positive constants. \end{example} Let $h=1/(N+1)$, where $N$ is an integer. Denote the numerical approximations of the steady state solution $u(t,x)$ at the interior points $x_j=jh$ by $u_j,j=1,2,\dots,N$. Using the central difference scheme $u_{xx} \approx (u_{j+1}+u_{j-1}-2u_j)/h^2$, the trapezoidal rule for the integral, and denoting $z=[u_1,\dots,u_N]^T$ and $y=v$ we obtain the following separable system for the steady state solutions of the differential equations: $$\begin{bmatrix} E(y,\mu) \\ B(y) \end{bmatrix} z + \begin{bmatrix} 0_{N \times 1} \\ (\mu-\mu_0)y-c_2y^2 \\ \end{bmatrix} = 0, \label{eq:4.3}$$ where $B(y)=c_3hy[1,\dots,1]$ and the $N \times N$ matrix $$E(y, \mu) = \begin{bmatrix} -2+h^2(\mu-c_1y) & 1 & {} & {} \\ 1 & -2+h^2(\mu-c_1y) & \ddots & {} \\ {} & \ddots & \ddots & 1 \\ {} & {} & 1 & -2+h^2(\mu-c_1y) \end{bmatrix}. \label{eq:4.4}$$ The matrix $A(y^*,\mu^*)$ has rank deficiency one when $(y^*,\mu^*)=(0,\mu_0)$ for $\mu_0=4(N+1)^2sin^2(\pi/(2N+2))$, as may be seen by the fact that $\operatorname{Null}(E(0,\mu_0))$ has dimension one. The point $(y^*,\mu^*,z^*)=(0,\mu_0,0)$ is a bifurcation point. Equation \eqref{eq:4.3} for $(y,\mu,z) \in \mathbb{R}^{N+2}$ reduces to a system of the form \eqref{eq:2.11} in $(y,\mu,w) \in \mathbb{R}^3$. As a numerical example we choose $N=19$, so $h=0.05$ and $\mu_0=9.84932752388982$. We choose $c_1=c_2=c_3=1$. Suppose we know a point $(\bar y, \bar \mu, \bar z) =(0.001, 9.848, (0.0002, \dots, 0.0002))$ near the bifurcation point $(y^*,\mu^*,z^*)=(0,\mu_0,0)$. The 19 singular values of $A(\bar y, \bar \mu)$ are $3.9507,3.8774,\dots,0.5612,0.3573,0.1934,0.0733,0.0002.$ The last number drops abruptly, which suggests that $A(y^*,\mu^*)$ has rank deficiency one. We obtain $\bar w =8.0361 \times 10^{-4}$ from $\bar{w}= \mathit{l}^T \bar{z}$. The singular value decomposition of $f'(\bar y, \bar \mu, \bar w)$ gives two singular values $\widetilde{\sigma}_1=0.0033$, $\widetilde{\sigma}_2=0.0000$ which implies $d=2$. Thus $\lambda^{(0)}=0 \in \mathbb{R}^2$. We randomly choose $\gamma=(0.3721, 0.6843)^T$. Table \ref{tb:2} gives data for three iterations and shows quadratic convergence. The computed approximation $z^*$ satisfies $\|z^*\|_2 \approx 2.2462 \times 10^{-14}$. If the underestimate $d=1$ is used then the rate of convergence is reduced to linear as mentioned above. \begin{table}[thbp] \caption{Newton iterations for computing a bifurcation point in Example \protect\ref{ex:2}} \label{tb:2} \begin{center} \begin{tabular}{||c||c|c|c|c||c||} \hline\hline $k$ & $y^{(k)}$ & $\mu^{(k)}$ & $w^{(k)}$ & $\lambda_1^{(k)}$ & $\|d^{(k)}\|_2$ \\ \hline\hline 0 & ~1.0000e-03 & 9.8480e+00 & ~8.0361e-04 & ~0.0000e+00 & 1.8461e-03 \\ 1 & -1.1537e-10 & 9.8493e+00 & ~1.3079e-09 & ~3.2664e-09 & 2.1661e-06 \\ 2 & -2.0761e-16 & 9.8493e+00 & -4.5227e-14 & -9.2603e-22 & 2.3200e-14 \\ 3 & -1.0311e-16 & 9.8493e+00 & -2.2462e-14 & ~2.8124e-30 & ----------- \\ \hline\hline \end{tabular} \end{center} \end{table} \begin{thebibliography}{99} \bibitem{beyn} {W.-J. Beyn}; \emph{Numerical methods for dynamical systems}, in Advances in Numerical Analysis, Vol. 1: Nonlinear Partial Differential Equations and Dynamical Systems, W. Light ed., Clarendon Press, Oxford, 1991, pp.~175--236. \bibitem{chow:shen} {S.-N. Chow and Y.-Q. Shen}; \emph{Bifurcations via singular value decompositions}, Appl. Math. Comp. 28(3) (1988), pp.~231--245. \bibitem{golu:leve} {G.~H. Golub and R.~J. LeVeque}; \emph{Extensions and uses of the variable projection algorithm for solving nonlinear least squares problems}, in Proceedings of the Army Numerical Analysis and Computing Conference, US Army Research Office, Washington DC, 1979, pp.~1--12. \bibitem{golu:pere:1} {G.~H. Golub and V. Pereyra}; \emph{The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate}, SIAM J. Numer. Anal. 10(2) (1973), pp.~413--432. \bibitem{golu:pere:2} {G.~H. Golub and V. Pereyra}; \emph{Separable nonlinear least squares: the variable projection method and its applications}, Inverse Problems 19 (2003), Topic Review, R1--R26. \bibitem{golu:vanl} {G.~H. Golub and C.~F. Van Loan}; \emph{Matrix Computations}, 3rd ed., Johns Hopkins, Baltimore, 1996. \bibitem{grie:redd} {A. Griewank and G.~W. Reddien}; \emph{Characterization and computation of generalized turning points}, SIAM J. Numer. Anal. 21 (1984), pp.~176--185. \bibitem{luke} {G.~G. Lukeman}; \emph{Separable Overdetermined Nonlinear Systems: An Applications of the Shen-Ypma Algorithm}, VDM Verlag Dr. M\"{u}ller, Saarbr\"{u}cken, 2009. \bibitem{morr} {D.~B. Morrison}; \emph{Application of the Shen-Ypma algorithm for nonlinear systems with some linear unknowns}, MSc Thesis, Dalhousie Univ., Halifax, Nova Scotia, 1998. \bibitem{mull:stok} {K. M. Mullen and I. M. van Stokkum}; \emph{The variable projection algorithm in time-resolved spectroscopy, microscopy and mass spectrometry applications}, Numer. Algorithm, 51(3)(2009), pp.~319-340. \bibitem{niko:bozh:nede:zlat} {S. Nikolov, B. Bozhov, V. Nedev and V. Zlatanov}; \emph{The Sherman system: bifurcation, regular and chaotic behaviour}, C. R. Acad. Bulgare Sci. 56(5) (2003), pp.~19-24. \bibitem{orte:rhei} {J.~M. Ortega and W.~C. Rheinboldt}; \emph{ Iterative Solution of Nonlinear Equations in Several Variables}, Academic Press, New York, 1970. \bibitem{rabi:redd} {P.~J. Rabier and G.~W. Reddien}; \emph{ Characterization and computation of singular points with maximum rank deficiency}, SIAM J. Numer. Anal. 23 (1986), pp.~1040--1051. \bibitem{shen} {Y.-Q. Shen}; \emph{Computation of a simple bifurcation point using one singular value decomposition nearby}, Computing 58(4) (1997), pp.~335--350. \bibitem{shen:2} {Y.-Q. Shen}; \emph{Computation of a Hopf bifurcation point using one singular value decomposition nearby}, in: Differential Equations and Applications, P. W. Bates, S.-N. Chow, K. Lu and X. Pan eds., International Press, Cambridge, 1997, pp.~277--285. \bibitem{shen:3} {Y.-Q. Shen}; \emph{Computation of a multiple bifurcation point using one singular value decomposition nearby}, Dynam. Cont. Disc. Impul. Syst. 6(1) (1999), pp.~53--68. \bibitem{shen:ypma} {Y.-Q. Shen and T.~J. Ypma}; \emph{Solving nonlinear systems of equations with only one nonlinear variable}, J. Comput. Appl. Math. 30(2) (1990), pp.~235--246. \bibitem{shen:ypma:2} {Y.-Q. Shen and T.~J. Ypma}; \emph{Solving separable nonlinear equations with Jacobians of rank deficiency one}, in LNCS 3314: Computational and Informational Sciences, J. Zhang, J.~-H. He and Y. Fu eds., Springer, New York, 2004, pp.~99--104. \bibitem{shen:ypma:3} {Y.-Q. Shen and T.~J. Ypma}; \emph{A uniform approach to computing dynamical equilibria}, Can. Appl. Math. Q. 14(3) (2006), pp.~343--359. \bibitem{shen:ypma:4} {Y.-Q. Shen and T.~J. Ypma}; \emph{Solving rank-deficient separable nonlinear equations}, Appl. Numer. Math. 57(5-7) (2007), pp.~609-615. \bibitem{shen:ypma:5} {Y.-Q. Shen and T.~J. Ypma}; \emph{Numerical bifurcation of separable parameterized equations}, Elect. Trans. Numer. Anal., 34 (2009), pp.~31-43. \bibitem{sher} {S. Sherman}; \emph{A third-order nonlinear system from a nuclear spin generator}, Contrib. Diff. Equation 2 (1963), pp.~197-227. \bibitem{watk} {D.~S. Watkins}; \emph{Fundamentals of Matrix Computations}, Wiley, New York, 1991. \bibitem{ypma:shen} {T.~J. Ypma and Y.-Q. Shen}; \emph{Solving N+m nonlinear equations with only m nonlinear variables}, Computing 44(3) (1990), pp.~259--271. \end{thebibliography} \end{document}