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

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

T = 1;


ux = 0;

u_num = zeros(n_t, n_x);
lamb = ones(n_t, n_x);

for i = 1:n_x
    
   u_num(n_t, i) = x_vec(i) -1; 
    
end


for count = 1:20
    
    for k = 1:n_t -1
        t = n_t -k;
        
        for i = 1:n_x
            
            % determin u
            
            
            
            if (lamb(t,i) < 0) || (i == 1)
                
               ux = ( u_num(t+1,i+1) - u_num(t+1,i) ) / Delta_x;
                
            else
            
            ux =  ( u_num(t+1,i) - u_num(t+1,i-1) ) / Delta_x;
                
            end
            
            
            u_num(t, i) = - sin(lamb(t,i)) * ux * Delta_t + u_num(t+1,i);
            
        
       % determine new lamb
       
       if ux > 0 
           lamb(t,i) = 0;
           
       else 
           
           lamb(t,i) = pi;
       
       end
        
    end
    
    end
    
    
end
 
 % Plot
% numerical solution
subplot(2,1,1);
zlim([-2,2]);
surf(t_vec,x_vec,u_num');

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