I coded a fully implicit finite difference model in MATLAB. It is helpful to discritize the vertical domain with uniform internode distances — rather than a domain with changing distances as deformation occurs — by rewriting the final form of the equation in the material coordinate, m, rather than in the Cartesian coordinate, z. That equation is not shown, but the steps you’d need to take are to (1) solve for dm/dz by differentiating the definition of the material coordinate, and (2) writing all d/dz terms as (dm/dz)*(d/dm) through the product rule.