For solving nonlinear partial equations, say $-\Delta y + f(y) = 0$, there are two common approaches:
Fixed point iteration: Pick $y^0$ and for $k=1,\dots$ solve
$$ -\Delta y^{k+1} + f(y^k) = 0.$$
(That is, you replace $y$ by $y^{k+1}$ in all linear terms and by $y^k$ in all nonlinear terms.)
Newton's method: Pick $y^0$ and for $k=1,\dots$ solve
$$ -\Delta \delta y + f'(y^k)\delta y = -(\Delta y^k+f(y^k)) \tag{1} $$
and set $y^{k+1} = y^k + \delta y$.
Which to use depends on the nonlinear term $f$; in the context of Reynolds equation, $f$ is quadratic, and this means that Newton's method is a good choice.
Now to the main question: How to actually compute $\delta y$. The easiest (and, at least for me, least error-prone way) is to apply your finite-difference discretization to $(1)$ and solve the discretized problem like you would any linear PDE.
Alternatively, let's say you already have a discretization of the nonlinear problem and want to solve the resulting nonlinear finite-dimensional system. From your question, I take it that you have a two-dimensional domain which you have discretized using $N:=mn$ grid points, resulting in matrices $K_h\in\mathbb{R}^{N\times N}$ (representing the linear part) and $F_h(y_h)\in\mathbb{R}^{N\times N}$ (representing the nonlinear part), where $y_h\in\mathbb{R}^N$ is the discrete solution. Since the dimensions only depend on the product $N$, this means that Newton's method does not care whether $n>m$ or not (your discretization will, of course).
Applying Newton's method to $K_h y_h + F_h(y_h) = 0$, you have to solve for $\delta y_h \in \mathbb{R}^N$ in
$$(K_h + \nabla F_h(y_h^k)) \delta y_h = - K_h y_h^k - F_h(y_h^k),$$
where $\nabla F_h(y_h)\in\mathbb{R}^{N\times N}$ is the Jacobi matrix of $F_h$ with entries $(\nabla F_h(y))_{ij} = \frac{\partial F_h(y_h)_i}{\partial (y_h)_j}$. This is a standard linear system, which you can solve with any of the usual methods. There is one additional point: If you pick an iterative method, you can exploit the fact that for $y^k$ far from the solution $y$, your linearization error will be large and likely dominate the approximation error, and thus there is no need to solve the linear system very accurately in the beginning. (Of course, you then need to increase the accuracy as you get closer to $y$.) This is known as an inexact Newton(-Krylov) method.
