M=4;
lambda=1/6;
N=M^2/2/lambda;
L=1;
T=0.5;
dx=L/M;
dt=T/N;
r=lambda;
for i=1:M+1
   v0(i)=(i-1)*dx+1000*sin(pi*(i-1)*dx);
end
for j=1:2
for i=2:M
   v1(i)=r*v0(i+1)+(1-2*r)*v0(i)+r*v0(i-1);
end
v1(1)=0;
v1(M+1)=1;
v0=v1;
end