 
delta_x = 0.2;
x = -1:0.2:1;
n_x = length(x);

delta_t = 0.2;
t_vec = 0:0.2:1;
n_t = length(t_vec);
T = 1;
H_1 = 0;
H_2 = 0;

v_num = rand(n_t, n_x);
alpha = ones(n_t,n_x);

for count = 1:20
    
    for s=1:n_t-1
    t = n_t - s ;
    for i = 2:n_x -1
        v_num
        H_1 = 1*(v_num(t,i) - v_num(t,i-1))/delta_x  -1;     % here is where alpha is equal to 1
        H_2 = (-1)*(v_num(t,i+1) - v_num(t,i))/delta_x  -1; % here is where alpha is equal to -1
        
        alpha(t,i) = 1;
        if H_2 > H_1
            alpha(t,i) = -1;
        end
        
        if alpha(t,i) == 1
            v_num(t,i) = (-v_num(t+1,i)/delta_t - v_num(t,i-1)/delta_x -1 )/(1/delta_t  -1/delta_x );
        else
             v_num(t,i) = (-v_num(t+1,i)/delta_t - v_num(t,i+1)/delta_x -1 )/(1/delta_t  -1/delta_x );
        end
        
        
    end
    
    end
    
    
end
 
 % Plot
% numerical solution
subplot(2,1,1)
surf(t_vec,x,v_num');

% exact solution
subplot(2,1,2)
v_sol = @(t,x) min(1-t, 1- abs(x) );
xmin = -1;
xmax = 1;
tmin = 0;
tmax = 1;
ezsurf(v_sol,[tmin,tmax,xmin,xmax]);