Edit: I'm going to give some more background and derivation to show how I got to these equations. I am basically following the derivation that is found in the appendix of the following paper:
R.M. Wood and B.E. Watts, "The Flow, Heat, and Mass Transfer Characteristics of Liquid Films on Rotating Disks", Trans. Instn Chem. Engrs., Vol 51, 1973.
The derivation is for the laminar fluid flow on the top of a rotating disk, where the liquid is dispensed on the center only.
Parameters: Fluid flow in the $r$-direction is $u(r,z)$, flow in the $\theta$-direction is $v(r,z)$. $\omega$ is the rotation rate of the disk, $h(r)$ is the height of the fluid film. Liquid flow rate is $Q=\int_0^{2\pi}\int_0^{h(r)}u(r,z)\,dz\,r\,d\theta=2\pi\,r\,h\,u_m$, where $u_m$ is the mean velocity in the $r$-direction, or $u_m(r)=\frac{1}{h(r)}\int_0^{h(r)}u(r,z)\,dz$. Similarly $v_m(r)=\frac{1}{h(r)}\int_0^{h(r)}v(r,z)\,dz$.
Assumptions: There is rotational symmetry such that $\frac{\partial}{\partial \theta}=0$. Flow is only horizontal so that flow in the $z$-direction $w=0$. Similarly $\frac{\partial P}{\partial z}=0$.
Assume that flow in the $r$-direction is parabolic and that $u=0$ at $z=0$ and $\frac{\partial u}{\partial z}=0$ at $z=h$. This gives us $u$ in terms of the local mean $u_m$ and fluid height $h$:
$$u(r,z)=3u_m\left[ \frac{z}{h}-\frac{1}{2}\left(\frac{z}{h} \right)^2 \right]$$
We do a similar assumption for flow in the $\theta$-direction, so that in this case $v=\omega\,r$ at $z=0$ and $\frac{\partial v}{\partial z}=0$ at $z=h$. This gives us:
$$v(r,z)=r\,\omega+\frac{3}{2}(r\,\omega- v_m)\left[\left(\frac{z}{h} \right)^2 -2 \frac{z}{h}\right]$$
Derivation
At this point we take the steady-state Navier-Stokes equations in cylindrical coordinates, assuming that $\frac{\partial P}{\partial r}=0$, keeping the advection terms for $u$ and $v$, and keeping the $\frac{\partial^2}{\partial z^2}$ viscous terms:
$$u\frac{\partial u}{\partial r}-\frac{v^2}{r}=\nu \left(\frac{\partial^2 u}{\partial z^2} \right)$$
$$u\frac{\partial v}{\partial r}+\frac{vu}{r}=\nu \left(\frac{\partial^2 v}{\partial z^2} \right)$$
At this point you integrate both equations from $0$ to $h$ in the $z$-direction. Using the above assumptions for $u(r,z)$ and $v(r,z)$, you can perform the integration for the viscous term on the RHS of the N-S equations. For the advection terms on the LHS, Wood and Watts then simply replace $u$ with $u_m$ and $v$ with $v_m$ and throw in an $O(1)$ constant for each, eliminating the need to solve the integral for the advection terms. Finally you use the definition of $Q=2\pi\,r\,h\,u_m$ to eliminate $h$ from the equations with $h=Q/(2\pi\,r\,u_m)$. For Wood and Watts' derivation this gives the following ODE's:
$$u_m \frac{du_m}{dr}=\frac{{v_m}^2}{r}-\frac{12 K_1 \nu\, \pi^2\,r^2\,{u_m}^3}{Q^2}$$
$$u_m \frac{dv_m}{dr}=-\frac{u_m v_m}{r}+\frac{12 K_2 \nu\, \pi^2\,r^2\,(r\,\omega-v_m){u_m}^2}{Q^2}$$
My goal however, is to eliminate the need for the $K_1$ and $K_2$ fudge factor constants by using the full assumed forms of $u(r,z)$ and $v(r,z)$ and integrating all the terms on the LHS of the NS equations. Finishing the integration and simplification (done using MATLAB symbolic math package) I have the following two equations:
$$69ru_m\frac{du_m}{dr}=8r^2\omega^2-16r\omega v-21{u_m}^2+48{v_m}^2-\frac{480\pi^2\nu r^3{u_m}^3}{Q^2}$$ $$48ru_m\frac{dv_m}{dr}=-21r(v_m-r\omega)\frac{du_m}{dr}-69u_m v_m+37r\omega u_m+\frac{480\pi^2\nu r^3(v_m-r\omega){u_m}^2}{Q^2}$$
Problem
I want to solve them simultaneously for $u(r)$ and $v(r)$, and initial conditions of $u(r_0)$ and $v(r_0)$ are known.
How to use the Runge-Kutta method is not my question. I've used it in the past and know how it works. My question/problem comes from the $\frac{du}{dr}$ term in the 2nd equation. Without that term I could just solve the two equations simultaneously using Runge-Kutta quite easily. But with it I'm having trouble understanding the best way to solve it. I could simply solve the first equation for $\frac{du}{dr}$ and substitute it into the second equation, but that would make the formulation really messy and I'd like to avoid that. Is there a way I could introduce a simple 3rd equation to deal with that? I'm thinking of something similar to how when you solve a 2nd degree ODE like $y''=f(x)$ using the Rung-Kutta method you split it into two equations:
$$y_1'=y_2$$
$$y_2'=f(t)$$
However I'm having trouble coming up with a formulation that would let me do that.
Update
With some help from the folks over at Computational Science Stackexchange, I have solved the system of equations. The solution can be seen here.
