function burgers ( )
clear

L = input('Area length : ');
U0 = input('Velocity amplitude : ');
nu = input('Viscosity : ');
Re = U0 * L / nu M = input('Numbrer of elements : ');
Tf = input('Final time : ');
dt = input('Time step : ');
dts = input('Backup time step : ');
dx = L / M; dxa = dx / L;
dta = dt * U0 / L;
dtsa = dts * U0 / L;
n = floor(Tf / dt);
ns = floor(Tf / dtsa); A = zeros(M + 1);
V = zeros(M + 1, ns + 1);
X = zeros(M + 1, 1);
U = zeros(M + 1, 1);
d = zeros(M + 1, 1); for i = 1 : M + 1
     X(i, 1) = (i - 1) * dx;
end for i = 1 : M + 1
     U( i, 1) = sin((i - 1) * 2 * pi / M);
     V( i, 1) = sin((i - 1) * 2 * pi / M);
end for j = 1 : n       for i = 2 : M
           A(i, i - 1) = - U(i - 1, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
           A(i, i) = 1 / dta + 1 / (Re * dxa * dxa);
           A(i, i + 1) = U(i + 1, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
           d(i, 1) = (U(i - 1, 1)^2 - U(i + 1, 1)^2) / (4 * dxa) + (U(i + 1, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(i + 1, 1) + (1 / dta - 1 / (Re * dxa * dxa)) * U(i, 1) + (- U(i - 1, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(i - 1, 1);
     end      A(1, M + 1) = - U(M + 1, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
     A(1, 1) = 1 / dta + 1 / (Re * dxa * dxa);
     A(1, 2) = U(2, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
     d(1, 1) = (U(M + 1, 1)^2 - V(2, 1)^2) / (4 * dxa) + (U(2, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(2, 1) + (1 / dta - 1 / (Re * dxa * dxa)) * U(1, 1) + (- U(M + 1, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(M + 1, 1);      A(M + 1, M) = - U(M, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
     A(M + 1, M + 1) = 1 / dta + 1 / (Re * dxa * dxa);
     A(M + 1, 1) = U(1, 1) / (4 * dxa) - 1 / (Re * 2 * dxa * dxa);
     d(M + 1, 1) = (U(M, 1)^2 - V(1, 1)^2) / (4 * dxa) + (U(1, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(1, 1) + (1 / dta - 1 / (Re * dxa * dxa)) * U(M + 1, 1) + (- U(M, 1) / (4 * dxa) + 1 / (Re * 2 * dxa * dxa)) * U(M, 1);      InvA = inv(A);
     U = InvA * d;      if (mod(j * dta, dtsa) == 0)
         for i = 1 : M + 1
               V(i, floor(j * dta / dtsa) + 1) = U(i, 1);
         end
     end
end V = V * U0
plot(X, V, 'r-');