Ah one more thing.first, you can evaluate the nonlinear terms in the current time step, solve the system with regular gauss elimination. Then use those values as the initial guess for the solution of the original equations where the nonlinear terms are evald in the forward time.
It seems fine to me. one issue with newton's method is you should be close enough to the root, otherwise it might not converge. maybe you might change the IC, if its ok, to something that satisfies the BC's also. then the guess values(the prev values) might be closer to the solution.
Can you show your equations you used for the boundary nodes?Also for inner nodes?i can check them for you.
How i do it is, use central diff to discretize bc's. You solve for the imaginary nodes (for 0, its one right outside say 0-, for 1 its 1+). Then you write the discrete form of the PDE then...