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=2: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=2:J+1+N-n
        u1(j)=u0(j)-(1/2)*(tau/h)*a*(u0(j+1)-u0(j-1))+(1/2)*(a^2)*((tau^2)/(h^2))*(u0(j+1)-2*u0(j)+u0(j-1));
    end
    for j=2:J+1+N-n
        u0(j)=u1(j);
    end
end
for j=2:J+1;
    axis_x(j)=x(j);
    appr_u(j)=u1(j);
end
plot(axis_x,appr_u)