Linearising Nonlinear Coupled Partial Differential Equations - Alfvénic Diffusion

I am trying to solve the following coupled partial differential equations with a finite difference scheme:
$$partial_tf+vpartial_zf+partial_zfrac{1}{W}partial_zf=0$$
$$partial_tW+vpartial_zW-partial_zf=0,$$

where both f and W depend on t and z.
In order to solve both equations simultaneously I am trying to use vector notation.
However the diffusion term in the first equation is non-linear.
How can I linearise this?

Background:

I have previously solved them in a semi-implicit scheme one after another, using the result of one equation in the other one. This however is explicit and requires very small timesteps, making the program to slow to be useful.

The basic idea is now to write them in a vector/matrix form and solve them simultaneously – again using finite differences.
This requires the equation to be of the form:
$$begin{array}{c}
f
W
end{array}=vec{f}$$

$$partial_tvec{f}=textbf{C}_0vec{f}+textbf{C}_zpartial_zvec{f}+textbf{C}_{zz}partial_{z}^2vec{f},$$

where all matrices $textbf{C}$ are independ of $vec{f}$.

The full problem consists of an additional dimension treated in LOD and a discontinuity in $v$ and $W$. Therefore, I do not expect standard ODE solvers to work.

Is there any way (substitution, or other) to make this possible?