It sounds like you're trying to solve a system of PDEs using FEM where two (or more equations) are dependent on each other, right?
If that's the case you have a couple of options:
1. "Staggered Solves": The idea here is to solve the equations one at a time and pass the solution from one solve into the next. You might iterate on this process (called Picard iteration) until you reach convergence.
2. "Fully Coupled" analysis where you use a Newton method to solve all of the equations simultaneously.
Both of these are equally easy to do using MOOSE ( http://mooseframework.org ). MOOSE was originally created to solve (using a "Fully Coupled" solver) systems of PDE's. Lately we've made large advances on doing staggered solves with MOOSE: http://mooseframework.org/blog/picard-with-multiapps/