On a cursory look it seems you are solving a advection diffusion reaction equations. The first one seems to be conservation of mass while the second one is AD equation for variable phi. You can always discretize using finite differences and use upwind schemes for advection terms,central differences for diffusion and a second order implicit or a semi implicit scheme like crank Nicolson for transient term should work. Since we have cross terms and it looks nonlinear a stability analysis of linearized version can help you derive a CFL like condition. I have written a second order accurate axi symmetric code for some industrial heterogeneous catalytic systems and the above said method worked okay though the source terms were a tad simpler. My understanding of FE is that it has some issues in handling advection terms and hence nearly all commerical fluid codes are erstwhile of FD type and now of Finite volume which uses a conservative form of naiver stokes equations.