I wish to numerically solve the following PDE. Although there are some complete discussions for solving PDEs in tutorial/NDSolvePDE, there is no hint for the nonlinear case by discretization. Thus, I will be thankful to receive some helps on the following NPDE where x in [0,1], t in [0,2],

M=8, NN=8, m = M - 1; n = NN - 1; alpha = 5.; beta = 4.; c = 0.05; T = 2.; h = (1. - 0.)/M; k =T/NN. The problem is Subscript[u, t] + u Subscript[u, xx] = c Subscript[u, xx] with I.C. Subscript[u, x, 0] = (2. c beta Pi Sin[Pi Subscript[x, i]])/(alpha + beta Cos[Pi Subscript[x, i]]) and B.C. Subscript[u, 0, t] = 0; Subscript[u, 1, t] = 0;

I tried the backward finite difference for Subscript[u, t] and the central FD for the others. I wrote the following code, but I think there are some gaps in it. Because, the approximate solutions do not match the exact one (u[x_, t_] := (2 c beta Pi Exp[-c Pi^2 t] Sin[Pi x])/(alpha + beta Exp[-c Pi^2 t] Cos[Pi x]);). Note that Subscript[w,i,j] stands for the approximation in the grid point (x_i,t_j).

 M = 8; NN = 8;
m = M - 1; n = NN - 1; alpha = 5.; beta = 4.; c = 0.05; T = 2.; h = (1. - 0.)/M; k =T/NN;
(*Defining the Grid points*)
Table[Subscript[x, i] = 0 + i h, {i, 0, M}]; Table[Subscript[t, j] = 0 + j k, {j, 0, NN}];
(*Defining the Initial Conditions*)
For[i = 1, i <= m, i++, Subscript[w, i, 0] = (2. c beta Pi Sin[Pi Subscript[x, i]])/(alpha + beta Cos[Pi Subscript[x, i]])];
(*Defining the Boundary Conditions*)
For[j = 1, j <= n, j++, Subscript[w, 0, j] = 0]; For[j = 1, j <= n, j++, Subscript[w, 1, j] = 0];
(*Defining the nonlinear equations due to discretization*)
For[i = 1, i <= m, i++,{For[j = 1, j <= n, j++,f[i, j] = Subscript[w, i, j] + (k/(2 h)) Subscript[w, i, j] (Subscript[w, i + 1, j] - Subscript[w, i - 1, j]) - (c k/(h^2)) (Subscript[w, i + 1, j] - 2 Subscript[w, i, j] + Subscript[w, i - 1, j]) - Subscript[w,i, j - 1]]}];
F = Flatten[Table[f[i, j], {i, 1, m}, {j, 1, n}]]; Dimensions[F];
F // MatrixForm;
Vec = Flatten[Table[Subscript[w, i, j], {i, 2, M}, {j, 1, n}]];
(*Finding the solutions*)
Sol = Part[NSolve[F, Vec, Reals], 1]

Any suggestion is appreciated. In fact, what would be the final nonlinear system of equations resulting of discretization?