I am trying to come up with a finite difference method for the equation: $$ u_t=-e^{-x}u_x+e^{-x}u-Le^{-x}u_z+Mu_z\left[Ef''(u)u_{zz}-\epsilon^2u_{zzz}\right] + Mu\left[Ef'''(u)(u_z)^2+Ef''(u)u_{zz}-\epsilon^2u_{zzzz}\right] $$ on a simple rectangle given by $D=[0,L]\times[0,1]$. The boundary conditions associated with the system are linear. As such, I was able to approximate them at the corresponding boundary using centered difference operators and then rearrange to solve for the "ghost points" immediately outside $D$. For example: $$ u_z=0 \ \ \text{on } z=1 \implies u_{j,M+1}=u_{j,M-1} $$ Note I am discretizing $D$ using points $(x_j,z_k)$: $$ x_j = j(\Delta x) \ \text{for } j=0,...,N \qquad\&\qquad z_k = k(\Delta z) \ \text{for } k=0,...,M$$ My question is in how to deal with the actual PDE. In the past when I've come across nonlinear equations, it was possible to rewrite the equation along the lines of something like: $$ u_t = L(u) + N(u) $$ representing the linear and nonlinear terms respectively. Then, in the past, I used an implicit method on the linear terms and then an explicit method on the nonlinear terms which worked quite well because in that case the highest order derivative was a linear term and so this greatly reduced the time-stepping criteria for the explicit method. This is clearly not the case here though.
Additionally, the function $f(u)$ is a cuartic polynomial and so I am a little concerned as to how to handle these terms. Will it be sufficient to approximate that at the previous timestep, i.e. $f(u)\approx f(u^n)$?
My thought is that this probably needs to be solved as a system of nonlinear algebraic equations at each time step. However, I've never quite done this before and am looking for some advice on the best way to go about this. My hope is to solve this using Matlab - I have used the built-in fsolve function to solve small systems, but am unsure the viability of this for a system that could be as large as $(128)(64)=8192$ nonlinear equations?