L = 10;
D = 0.5;
s = 0.03;
Tmax = 2;
xdim = 25; tdim = 75;
x = linspace(0,L,xdim);
t = linspace(0,Tmax,tdim);
dx = x(2)-x(1); dt = t(2)-t(1);
q = dt/dx^2;
r1 = 0.75*L; r2 = 0.8*L;
u0 = zeros(1,xdim);
u0(find(x>=r1 & x<=r2)) = 1;
xDat = 2:xdim-1;
tDat = tdim;
nxDat = length(xDat);
ntDat = length(tDat);
Z = heat(D,u0,q,tdim);
u = Z(tDat,xDat);
uDat = u + s*randn(ntDat,nxDat);
figure(1); surf(x,t,Z); hold on;
if ntDat>1, mesh(x(xDat),t(tDat),uDat);
else set(plot3(x(xDat),t(tDat)*ones(1,nxDat),uDat,'r-o'),'LineWidth',3);
end; hold off; drawnow
N = 10000; m = 100; XD = 1; X = zeros(1,N); X(1) = XD;		
Z = heat(XD,u0,q,tdim); 
u = Z(tDat,xDat);
oLLkd = sum(sum(-(u-uDat).^2))/(2*s^2); 
LL = zeros(1,N); LL(1) = oLLkd;
w = 0.1;
for n = 2:N
 XDp = XD + w*(2*rand-1);
 if XDp > 0
   Z = heat(XDp,u0,q,tdim);      
   u = Z(tDat,xDat);
   nLLkd = sum(sum(-(u-uDat).^2))/(2*s^2);       
   alpha = exp(nLLkd - oLLkd);
   if rand < alpha
     XD = XDp;
     oLLkd = nLLkd;
     CZ = Z;
   end
 end
 X(n) = XD; LL(n) = oLLkd;
 if rem(n,m)==0
   figure(2); plot(X(1:n)); drawnow;
   figure(3); surf(x,t,CZ); hold on;
   if ntDat>1, mesh(x(xDat),t(tDat),uDat);
   else set(plot3(x(xDat),t(tDat)*ones(1,nxDat),uDat,'r-o'),'Linewidth',3);
   end; hold off; drawnow
   disp([N/m,n/m]);
 end
end
figure(2); plot(X);