N  = 200;
N1 = N+1;
N2 = round(N/2);
h  = 1/N;
x  = linspace(0,1,N+1);
mu = 1.0e-3;
ustar = [zeros(N2,1);ones(N1-N2,1)];
v = ustar + 0.1*randn(N+1,1);

D = diag([1;2*ones(N-1,1);1]) ...
  - diag(ones(N,1),+1) ...
  - diag(ones(N,1),-1);
u = (mu/h^2*D + eye(N+1))\v;

plot(x,ustar,'b',x,v,'r',x,u,'g')
legend({'u*', 'v', 'u'})

