• The final equation can be written as :  where A (resp. d) is a square matrix (resp. a vector) that is a function of Un. Finally we obtain : .
• Matlab is a mathematical language able to easily reverse a matrix. There is no use to write a subscript computing A-1.
• The goal of this program is to plot U for several values of the time on the same graphe in order to follow the shock evolution as it is schown here.
• V is a matrix where U is saved for chosen values of the time. At the end of the program V versus X is plotted.

function burgers ( )
clear

• Physical variables input :
L = input('Area length : ');
U0 = input('Velocity amplitude : ');
nu = input('Viscosity : ');
Re = U0 * L / nu
• Numerical variables input :
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);
• Matrix and vector definition :
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);
• Abscisse vector calculation :
for i = 1 : M + 1
X(i, 1) = (i - 1) * dx;
end
• Initial condition
for i = 1 : M + 1
U( i, 1) = sin((i - 1) * 2 * pi / M);
V( i, 1) = sin((i - 1) * 2 * pi / M);
end
• Velocity computation at t = j  * dt :
for j = 1 : n
•  A and d calculation :
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
• Boundary conditions, x = 0 :
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);
• Boundary conditions, x = L :
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);
• Velocity computation :
InvA = inv(A);
U = InvA * d;
•  Velocity backup :
if (mod(j * dta, dtsa) == 0)
for i = 1 : M + 1
V(i, floor(j * dta / dtsa) + 1) = U(i, 1);
end
end
end
• U plotting :
V = V * U0
plot(X, V, 'r-');