n=10000; th(1)=0; pdf=inline('exp(-x^2/2)');
    for i=2:n             thstar = th(i-1) + 3*(1-2*rand);
           alpha = min([1,pdf(thstar)/pdf(th(i-1))]);
    if    rand <= alpha     th(i)=thstar;
    else  th(i)=th(i-1);  end 
    hist(th,100)
end
