function x=hastings(fn,z0,d,nsim)
% runs the hasting=meetropolis algorithm to generate a transition probability density given by the function fn. 
%z0 is the initial value.d=diameter of random movement, nsim=number of simulations
%e.g. run as x=hastings('wigglyfn',[0 0],1,10000); plot(x(:,1),x(:,2),'.')
fold=eval(strcat(fn,'(z0(1),z0(2))'));   x=[]; 
N=ceil(1.4*nsim);
R=(2*rand(2,N)-1);  v=(sum(R.^2)<=1); R=R(:,v);    %generates random points in circle of radius 1.
for i=1:nsim
z=z0+d*(R(:,i))';
fnew=eval(strcat(fn,'(z(1),z(2))'));
if (rand(1,1)<(fnew/fold))
    z0=z;
    fold=fnew;
    x=[x; z];
else 
    x=[x;z0];
end
end