I want to numerically solve for $y(x)$. I have the boundary conditions $y(x_a)=1$ and $y(x_b)=1$ as well as the following second-order set of differential equations:
$$ \frac{\text{d}z(x)}{\text{d} x} = - \frac{f_1(x)}{x^2 f_2(x)} y(x) $$ $$ z(x) = \frac{ \frac{\text{d}y}{\text{d}x} + f_3(x) }{ c x^2 f_1(x) } $$
I understand that to solve a BVP I need to reduce these equations to first-order. How can I do that?
Edit: or are they already first order, I just need to rearrange the $z(x)$ equation?
To solve a BVP you do not need to reduce it to first order. Since your boundary conditions are both on $y(x)$, its much easier to get a second order equation for $y(x)$ and solve that. From $y(x)$ you can easily find $z(x)$ so your problem will be solved.
The trick to get a nice form for solving your problem is to take the derivative of your second equation and put the first one into the second : $$ -\frac{f_1(x)}{x^2f_2(x)}y(x) = \frac{d}{dx}\frac{\frac{dy}{dx}+f_3(x)}{c x^2 f_1(x)} $$ This is now a second order equation for $y(x)$, which you can solve using the usual tricks. I would recommend using finite differences if you're writing your own solver. Consider solving $y(x)$ on grids where $\vec{x} = [x_0, x_1, ..., x_{i-1}, x_i, x_{i+1}, ..., x_n]$ and $\vec{y} = [y_0, y_1, ..., y_{i-1}, y_i, y_{i+1}, ..., y_n]$. The operator $y''(x) = h^{-2}/2 *(y_{i-1} -2 y_i + y_{i+1})$ where $h = \Delta x$, the grid step size, can be easily represented by a sparse tridiagonal matrix (and the same goes for other operators, you can find tables for finite differences on wiki). You can enforce your boundary conditions in a matrix by putting $1$ on the main diagonal and you will end up looking at a matrix system that looks like $M \vec{y} = \vec{b}$ which can be solved by linear algebra techniques.
EDIT : On non-uniform grids, you can still use finite differences, but you have to slightly change the formulas. Typically, finite difference schemes are obtained from Taylor expensions $f(x_i + \Delta x_i) = f(x_i) + f'(x_i) \Delta x_i + O(\Delta x_i^2)$ and if your grid is non-uniform, the coefficients will depend on the local grid size. You can find such formulas easily on the internet (see http://fundamentalthinking.blogspot.com/2014/08/second-order-finite-difference-schemes.html for instance).
The matrix $M$ in your case does not depend on $y$ as your system is linear in $y$. Consider solving the equation $y'' + a(x) y' + b(x) y =c(x)$, a general second order ODE. In a linearly spaced grid, $y'(x) = (2h)^{-1} (y(x+\Delta) - y(x-\Delta))$, which translates into $y'_i = (2h)^{-1} (y_{i+1} - y_{i-1})$. The matrix $A$ for the operator $a(x) y'(x)$ is a $N*N$ matrix, one line for each element of $\vec{y}$. The $i^{th}$ line elements are given by $A_{i,i-1} = -(2h)^{-1}a(x_i)$ and $A_{i, i+1} = (2h)^{-1} a(x_i)$. Also notice $y(x) = I$ and $c(x)$ is actually a vector. This matrix construction is fine for all elements with the exception of the end points where you will enforce boundary conditions.