I just knew how to do Newton-Raphson iteration in time-independent 1D nonlinear differential equation. Then I applied to time-dependent 1D nonlinear differential equation, and I got confused.Below is just test time-dependent 1D nonlinear differential equation I want to solve (I just made this up from heat equation, to make it nonlinear)$$\frac{\partial u}{\partial t}+\frac{\partial^2 u}{\partial x^2}+u^2=0$$For this equation, I assumed I would need initial condition, and boundary condition$$u(x,t=0)=\left\{\begin{aligned}1\qquad&(-1\leq x\leq 1)\\0\qquad&\text{otherwise}\end{aligned}\right.,\qquadu(-5,t)=0,\quad u(5,t)=0$$I set the $x$ range of $-5 \leq x \leq 5$, and planned to create code in Matlab.I started with space discretization by using second-order difference$$\frac{\partial u}{\partial t}+\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}+u_{i}^2=0$$where $i$ denotes space discretization numberThen I applied Backward-Euler discretization for time discretization, which is$$\frac{d q(x,t)}{dt}=f(t,x,q(x,t))\qquad\longrightarrow\qquad\frac{q_{j+1}-q_{j}}{\Delta t}=f_{j+1}$$where $j$ denotes time discretization number(this is just intermediate step for me. I am planning to apply Trapezoidal Rule - Second order Backward Difference Formula (TR-BDF2) later)Now applying this, the equation looks like$$\frac{u_{i,j+1}-u_{i,j}}{\Delta t}+\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}+u_{i,j+1}^2=0$$Again, $i$ and $ j$ denotes space and time discretization respectively.Then I applied Newton-Raphson(NR) scheme, which is for give equation$$F(u)=0$$the solution is determined by iteration$$F'(u^{k})\delta u^{k}=-F(u^{k}), \qquad \qquad u^{k+1}=u^{k}+\delta u^{k}$$(I'm running out of space for discretization) where k denotes NR iteration number, and if $F(u)$ was system of equations, $F'(u)\delta u$ should be Jacobian.Do I have to set initial and boundary condition for $\delta u_{i,j}$ as well?Now my equation looks like$$\frac{\delta u_{i,j+1}^{k}-\delta u_{i,j}^{k}}{\Delta t}+\frac{\delta u_{i+1,j+1}^{k}-2\delta u_{i,j+1}^{k}+\delta u_{i-1,j+1}^{k}}{\Delta x^2}+2u_{i,j+1}^{k}\delta u_{i,j+1}^{k}=-\frac{u_{i,j+1}^{k}-u_{i,j}^{k}}{\Delta t}-\frac{u_{i+1,j+1}^{k}-2u_{i,j+1}^{k}+u_{i-1,j+1}^{k}}{\Delta x^2}-(u_{i,j+1}^{k})^2$$This is where I got stuck. How do I proceed after this? Am I missing something? like another boundary/initial conditions?Because I only know $u_{i,1}^{0}$ and $\delta u_{i,1}^{0}$. There are like, 11 unknowns due to backward differences ($j+1$).Some literature says that I should solve this equation for each time stepOr, is this right way to apply Backward-Euler scheme and NR scheme to the time-dependent nonlinear differential equation?Is there any good example solving time dependent nonlinear differential equation with Newton-Raphson iteration?
Implementation of Backward-Euler scheme, Newton-Raphson iteration scheme to time dependent nonlinear differential equation
nonlinear equations;newton method
You seem to be confused with what are unknowns in your algebraic equations (as mentioned also in the comment by  origimbo).Let say you have $N$ inner nodes for space discretization, i.e. $\Delta x = 10/(N+1)$. The values $u_{i,1}$ shall be obtained from initial condition (a more usual notation is $u_{i,0}$).The notion you have to solve your algebraic equations in each time step means that in each time step the values $u_{i,j}$ are known (either from initial condition or from previous time step) and your unknowns (the values to be found) are only $u_{i,j+1}$. Now this is the task for NR method.When applying NR method, you have to choose the values $u^0_{i,j+1}$, typically you take $u^0_{i,j+1}=u_{i,j}$. It means $F(u)=0$ represents $N$ nonlinear algebraic equations having $N$ unknowns $u_{1,j+1}$ up to $u_{N,j+1}$ having the form$$\frac{\delta u_{i,j+1}^{k}}{\Delta t}+\frac{\delta u_{i+1,j+1}^{k}-2\delta u_{i,j+1}^{k}+\delta u_{i-1,j+1}^{k}}{\Delta x^2}+2u_{i,j+1}^{k}\delta u_{i,j+1}^{k}=-\frac{u_{i,j+1}^{k}-u_{i,j}}{\Delta t}-\frac{u_{i+1,j+1}^{k}-2u_{i,j+1}^{k}+u_{i-1,j+1}^{k}}{\Delta x^2}-(u_{i,j+1}^{k})^2$$In above, the values $u_{i,j}$, $u_{*,j+1}^k$ are known, the unkowns in above are only three (but in each equation different ones), namely $\delta u_{i-1,j+1}^k$, $\delta u_{i,j+1}^k$, and $\delta u_{i+1,j+1}^k$. The boundary conditions are $\delta u_{0,j+1}^k=0$ and $\delta u_{N+1,j+1}^k=0$.By the way, the sign in your PDE before the diffusion term is wrong, it should be minus.