 

x_delta = 1/10;
t_delta = 1/10;

x = -1:x_delta:1;
x_inner = x(2:end-1);
n_x = length(x);
t = 0:t_delta:1;
n_t = length(t);

v_num = zeros(n_t,n_x-2);


eps = x_delta/2;

E_mo = zeros(n_x-2, n_x-2); % mo stands for minus one
E_po = zeros(n_x-2, n_x-2); % po stands for plus one

% Determine matricies E^{\alpha} for \alpha =1 and -1
for l = 1:n_x-2
    for j = 1:n_x-2
        
        if j == l-1
            
            E_mo(l,j) = 1/(2*x_delta) - eps/(x_delta^2);
            E_po(l,j) = -1/(2*x_delta) - eps/(x_delta^2);
            
        elseif j == l
            
            E_mo(l,j) = (2*eps)/(x_delta^2);
            E_po(l,j) = (2*eps)/(x_delta^2);
            
        elseif j == l+1
            
            E_mo(l,j) = -1/(2*x_delta) - eps/(x_delta^2);
            E_po(l,j) =  1/(2*x_delta) - eps/(x_delta^2);
            
        end
        
    end
end


F_alph = 1; % Its always constant 1 in this example, doesnt depent on the control or time


v_num(n_t,1:end) = 0; % v(T,.)=0
E_mobay = zeros(n_x-2,1);
E_pobay = zeros(n_x-2,1);

for k_bay = 1:n_t-1
    s = n_t-k_bay;
    
    E_mobay(:,1) = E_mo*(v_num(s+1,:)');
    E_pobay(:,1) = E_po*(v_num(s+1,:)');
   
    for l = 1:n_x-2
        maxpo = E_pobay(l,1) -F_alph;
        maxmo = E_mobay(l,1) -F_alph;
        
        v_num(s,l) = -t_delta*max(maxpo,maxmo) + v_num(s+1,l);
        
    end
    
end

% Plot
% numerical solution
subplot(2,1,1)
surf(t,x_inner,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]);