Assume the equations are discretized on the $\tau$ grid, and introduce several column-vectors, using the notation where the superscript stands for the grid point index i $\in$ [1,...,n].
$\vec{\tau}=\left[ \tau^1,\tau^2,...,\tau^{n} \right] $
$\vec{\phi_1}=\left[ \phi_1^1,\phi_1^2,...,\phi_1^{n} \right] $
$\vec{\phi_2}=\left[ \phi_2^1,\phi_2^2,...,\phi_2^{n} \right] $
$\vec{f_{11}}=\left[ f_{11}^1,f_{11}^2,...,f_{11}^{n} \right] $
$\vec{f_{12}}=\left[ f_{12}^1,f_{12}^2,...,f_{12}^{n} \right] $
$\vec{f_{21}}=\left[ f_{21}^1,f_{21}^2,...,f_{21}^{n} \right] $
Next, the second derivative operator on the grid is a matrix
$\hat{\Delta}$ of size $n \times n$.
Also, use $\alpha (k) = c_2(k) / c_1(k)$, and denote
$c_1(k) = \lambda$ and $c_2(k) = \alpha \lambda$
Now, let's write both equations as a single linear system of size
$2n \times 2n$ for the compound state column-vector $\vec{\phi}$ which is
\begin{equation}
\vec{\phi} = [\phi_1^1,\phi_1^2,...,\phi_1^{n},\phi_2^1,\phi_2^2,...,\phi_2^{n}]
\end{equation}
We'll use here several matrices of size $2n \times 2n$:
The left hand side matrix $L$
\begin{bmatrix}
\hat{I} & \hat{0}\\
\hat{0} & \alpha \hat{I}
\end{bmatrix}
the differential operator matrix $D$
\begin{bmatrix}
\hat{\Delta} & \hat{0} \\
\hat{0} & \hat{\Delta}
\end{bmatrix}
the cross-coupling matrix $C$
\begin{bmatrix}
\hat{0} & \hat{I} \vec{f_{12}} \\
\hat{I} \vec{f_{21}} & \hat{0}
\end{bmatrix}
and the forcing term matrix $F$
\begin{bmatrix}
\hat{I} \vec{f_{11}} & \hat{0} \\
\hat{0} & \hat{0}
\end{bmatrix}
Here $\hat{I}$ is the unit matrix of size $n \times n$,
$\hat{0}$ is the null matrix of size $n \times n$.
Now our problem can be written as
\begin{equation}
\lambda L \vec{\phi} =
\left(
D + C + F
\right) \vec{\phi}
\end{equation}
or
\begin{equation}
\lambda \vec{\phi} =
L^{-1} \left(
D + C + F
\right) \vec{\phi}
\tag{*}\label{*}
\end{equation}
This a linear eigenvalue problem, where the right hand side matrix contains the parameter $\alpha$, so the eigenvalues are functions of $\alpha$, which in turn is a function of $k$.
For a given value of $k$ there is a spectrum of eigenvalues, $\lambda_0, \lambda_1$ etc. Let's say we are interested in the smallest eigenvalue $\lambda_0$ which we'd just call $\lambda$. Solution of the eigenvalue problem
$\eqref{*}$ by standard methods of linear algebra defines a function $\lambda(k)$, on the other hand by construction $\lambda$ must be equal to $c_1(k)$. The root of the equation $\lambda(k) = c_1(k)$ (if it exists) defines the solution of the problem. Since $\lambda(k)$ and $c_1(k)$ are two nonlinear functions the root of the equation $\lambda(k) = c_1(k)$ has to be sought by some iterative numerical technique.