Trying to model a substance settling in water using an advection equation?

Question

I am trying to model a substance dispersed in a container of water gradually settling at the bottom. I am considering only one dimension. The top is at $z = 1$, and the bottom is at $z = 0$.

So at $t = 0$ the concentration $u = 0.2$ everyhwere.

As time goes on the substance moves to the bottom where it settles and eventually the concentration should be

$ u(z, t)=\begin{cases} 0 & \text{if} \ z > 0.2\\ 1 & \text{if} \ z <= 0.2\\ \end{cases} $

I thought an advection type equation might work so I used

$$\frac{\partial u}{\partial t} = - v(u)\frac{\partial u}{\partial z}$$

with initial condition

$ u(z, 0)=\begin{cases} 0 & \text{if} \ z = 1\\ 0.2 & \text{if} \ 0 < z < 1\\ 1 & \text{if} \ z = 0 \\ \end{cases} $

But while the concentration in the upper part of the container does go to $0$, the concentration in the bottom stays at $0.2$ instead of going to $1$.

I used the method of lines with finite differences. I've attached the code below so that hopefully someone can see where I've gone wrong and tell me how to fix it so that the concentration ends being $1$ between $z = 0$ and $z = 0.2$?

function mol_edbl
clear
clc

global n

% Time and spatial steps
n = 51;

% Time start and finish
t0 = 0.0;
tf = 1;
t_grid_points = linspace(t0,tf,n);

% Spatial start and finish
z0 = 0;
zf = 1;
z_grid_points = linspace(z0,zf,n);

dz = (zf-z0)/(n-1);

% Initial condition
for i=1:n
    u_l0(i) = 0.2;
    v_l0(i) = 0;
end

u_l0(1) = 1;

u_l0(n) = 0;


ics = [ u_l0 v_l0 ];

reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t, solved]=ode15s(@solve, t_grid_points, ics, options);

u_l_solved = solved(1:n,1:n);
v_l_solved = solved(1:n,n+1:2*n);

% Store numerical and analytical solutions, errors at z = 1/2
  n2=(n-1)/2.0+1;
  sine=sin(pi/2.0*0.5);
  for i=1:n
    u_l_plot(i) = u_l_solved(i,n);
    v_l_plot(i) = v_l_solved(i,n);
  end
%

% --------------------%
% Plot
% --------------------%
% **************************************** %
% Display Results - Solution in 3d - Style 2
% **************************************** %
figure(1);
set(gcf,'Position',[200 100 1400 700]);
subplot(1,2,1)

mesh(t_grid_points, z_grid_points, u_l_solved');
title('u_l');
xlabel('t, time')
ylabel('z, distance')
zlabel('u_l(z,t)')

subplot(1,2,2)
mesh(t_grid_points, z_grid_points, v_l_solved');
title('v_l');
xlabel('t, time')
ylabel('z, distance')
zlabel('v_l(z,t)')
rotate3d on;

% --------------------%
% Solver
% --------------------%
function dt = solve(t, u)

u_l = u(1:n);
v_l   = u(n+1:2*n);

for i=1:n
    dv_l(i) = 0;

    if i == 1
        du_l(i) = 0;
    elseif i == n
        du_l(i) = 0;
    else
        du_l_dz = (u_l(i+1) - u_l(i))/dz;
        v = - 1 + 2 * u_l(i);
        du_l(i) =  - v * du_l_dz;
    end
end

dt = [ du_l' ; dv_l' ];
end

end

u_l0(1) = 1;

u_l0(n) = 0;


ics = [ u_l0 v_l0 ];

reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t, solved]=ode15s(@solve, t_grid_points, ics, options);

u_l_solved = solved(1:n,1:n);
v_l_solved = solved(1:n,n+1:2*n);

% Store numerical and analytical solutions, errors at z = 1/2
  n2=(n-1)/2.0+1;
  sine=sin(pi/2.0*0.5);
  for i=1:n
    u_l_plot(i) = u_l_solved(i,n);
    v_l_plot(i) = v_l_solved(i,n);
  end
%

% --------------------%
% Plot
% --------------------%
% **************************************** %
% Display Results - Solution in 3d - Style 2
% **************************************** %
figure(1);
set(gcf,'Position',[200 100 1400 700]);
subplot(1,2,1)

mesh(t_grid_points, z_grid_points, u_l_solved');
title('u_l');
xlabel('t, time')
ylabel('z, distance')
zlabel('u_l(z,t)')

subplot(1,2,2)
mesh(t_grid_points, z_grid_points, v_l_solved');
title('v_l');
xlabel('t, time')
ylabel('z, distance')
zlabel('v_l(z,t)')
rotate3d on;

% --------------------%
% Solver
% --------------------%
function dt = solve(t, u)

    u_l = u(1:n);
    v_l   = u(n+1:2*n);

    for i=1:n
        dv_l(i) = 0;

        if i == 1
            du_l(i) = 0;
        elseif i == n
            du_l(i) = 0;
        else
            du_l_dz = (u_l(i+1) - u_l(i))/dz;
            v = - 1 + 2 * u_l(i);
            du_l(i) =  - v * du_l_dz;
        end
    end

    dt = [ du_l' ; dv_l' ];
end

end