u=@(x) sin(pi*x);                                  % solution and corresponding
f=@(x) -pi^2*sin(pi*x);                            % right hand side
d=0.5;                                             % perturbation
NN=10*2.^(0:5);                                    % of the mesh, 0<=d<=0.5
err=[];
for j=1:length(NN)                                 % mesh refinement
  N=NN(j);
  if j==1
    h=1/N; xx=(0:h:1);
    xx(2:end-1)=xx(2:end-1)+(rand(1,N-1)-1/2)*d*h; % perturb primal mesh
    hhd=diff(xx);
    xxd=xx(1:end-1)+(1/2+(rand(1,N)-1/2)*d).*hhd;  % perturb dual mesh
    hh=diff([0 xxd 1]);
    figure(1);plot(xx,1,'o',xxd,1,'+');
  else
    xx(1:2:2*length(xx)-1)=xx;                     % refine mesh
    xx(2:2:end-1)=xxd;
    hhd=diff(xx);
    xxd=xx(1:end-1)+1/2.*hhd;
    hh=diff([0 xxd 1]);
  end;
  A=sparse(N-1);
  for i=1:N-1                                      % assemble matrix by going
    A(i,i)=-(1/hhd(i)+1/hhd(i+1))/hh(i+1);         % over all cell interfaces
    if i>1, A(i,i-1)=1/hhd(i)/hh(i+1); end;
    if i<N-1, A(i,i+1)=1/hhd(i+1)/hh(i+1); end;
  end;
  b=f(xx(2:end-1)'); ue=u(xx');
  ua=A\b;
  err(j)=max(abs(ue(2:end-1)-ua))
end
figure(2);loglog(1./NN,err,'--',1./NN,1./NN,'-',1./NN,1./NN.^2,'-')
