a=1.0;
h=0.1;
lambda=0.5;
tau=lambda*h;
T=0.5;
N=T/tau;
J=(1.0-0)/h+1;
% for initialization
for j=1:J+1+N
    x(j)=-(j-1)*h+1.0;
    if(x(j)>0.4&x(j)<0.6)
        u0(j)=1.0;
    else
        u0(j)=0;
    end
end
for n=1:N
    for j=1:J+1+N-n
         u1(j)=u0(j)-a*lambda*(u0(j)-u0(j+1));
    end
    for j=1:J+1+N-n
        u0(j)=u1(j);
    end
end
for j=1:J+1;
    axis_x(j)=x(j);
    appr_u(j)=u1(j);
end
plot(axis_x,appr_u)