I am trying to decide whether my derived PDE is linear - I am drawing a blank. If someone could demonstrate why it is/isn't that would be great! I have used the MHD equations of motion and induction in cylindrical polars to dervie the following. We assume that $\dfrac{\partial}{\partial\theta}=0$.
$\left(\dfrac{\partial}{\partial t}+\dfrac{\partial v_z}{\partial z}+v_z\dfrac{\partial}{\partial z}\right)W_\theta=\dfrac{1}{\mu}\dfrac{\partial}{\partial z}(B_\theta B_z v_\theta)-\dfrac{B_\theta^2}{2\mu}\dfrac{\partial v_z}{\partial z}$
I have that $W_\theta=\dfrac{v_\theta^2 \rho}{2}+\dfrac{B_\theta^2}{2}$ is the energy density, B is the magnetic field in three dimensions v is the velocity in three dimensions. All of the components only comprise of Z as I have used the definition for the area of a coronal hole ($A(z)=az^n$ for n=2) - this allows us to represent the plasma (solar wind) moving with a coronal hole.
This equation seems to be a particular case of the Cauchy momentum equation: $$ \rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} \right)= \nabla \cdot \mathbf{\sigma} + \mathbf{g} $$ in which $\sigma$ is the stress tensor (involving pressure and viscosity) and $\bf{g}$ is the body force (gravitation, magnetism, etc.). In this form, even in the absence of any term in RHS, the equation is non-linear, due to the $\bf{v} \cdot \nabla \bf{v}$.
It don't make that sense to talk about independent variables without proper context; in three dimensions (the most general case) this (vectorial) equation are actually three (scalar) equations, and there is four (or even five) variables: the three components of $\bf{v}$ and the pressure, that is hidden inside $\sigma$ (and, eventually, $\rho$, if your flow is compressible). You need also a mass conservation equation and a equation of state (and even an energy conservation equation, if this process is relevant).
Another source of non-linearity can be inside $\sigma$ if the viscous stress tensor is not Newtonian. In the most general case, each component of a Newtonian stress tensor is a linear combination of gradients of the velocity components. Any other viscous behavior is necessarily non-linear.
Returning, finally, to your case, the sources of non-linearity are the convective term $v_z \frac{\partial W_\theta}{\partial z}$ and the definition of $W_\theta$ itself. The first one is the same intrinsic non-linearity of Cauchy momentum equation and of the well-known Navier Stokes equation.