Numerical method needed for the following system of ODEs

tchar

I am required to solve the following system of ODEs numerically. Could someone suggest an appropriate methodology. These equations are phenomenological equations derived from irreversible thermodynamics. I have to solve for the flux terms given on the L.H.S. The coefficients of differentials on the R.H.S are constants and are known.


J_q= -L_11/T dT/dx-L_12 (dμ_HCL)/dx-L_13 ((d∅)/dx)
J_HCL= -L_21/T dT/dx-L_22 (dμ_HCL)/dx-L_23 ((d∅)/dx)
j= -L_31/T dT/dx-L_32 (dμ_HCL)/dx-L_33 ((d∅)/dx)