# Set up a tridiagonal for a system of equations

I am trying to solve the system of simultaneous equations:
$$\frac{\partial \rho_f\phi}{\partial t}+\frac{\partial}{\partial z}(\rho_f\phi v_f)= \frac{\partial F}{\partial t}$$
$$\frac{\partial \rho_s(1-\phi)}{\partial t}+\frac{\partial}{\partial z}(\rho_s(1-\phi) v_s)=-\frac{\partial F}{\partial t}$$
and
$$v_f-v_s=-\frac{k}{\phi \mu}\frac{\partial P}{\partial z}$$
where v_s and v_f are unknown except v_s is known at the lower boundary in z.
I will use explicit finite differences such that, I think, the discretized forms are
$$\rho_f^n\phi^n-\rho_f^0\phi^0+\frac{\Delta t}{\Delta z}[\rho_f^n\phi^nv_f^n-\rho_f^0\phi^0v_f^0]=F^n-F^0$$
$$\rho_s^n(1-\phi^n)-\rho_s^0(1-\phi^0)+\frac{\Delta t}{\Delta z}[\rho_s^n(1-\phi^n)v_s^n-\rho_s^0(1-\phi^0)v_s^0]=-(F^n-F^0)$$
and
$$v_f-v_s=-\frac{k}{\phi \mu}\frac{P_{z+1}-P_z}{\Delta z}$$
where superscripts n and 0 are for new time and current time, respectively.
But I do not understand what the strategy should be for building a tridiagonal matrix with these equations. What should I be doing? How can I set this up?

More on terms: v_f and v_s are velocities, rho_f and rho_s are densities, \phi is volume of material moving at velocity v_f, F is a mass flux, and k and
\mu are constants.