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 initial condition

Subscript[u, x, 0] = (2. c beta Pi Sin[Pi Subscript[x, i]])
                   / (alpha + beta Cos[Pi Subscript[x, i]]) 

and boundary condition

Subscript[u, 0, t] = 0; Subscript[u, 1, t] = 0;

I tried the backward finite difference (FD) 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?