function u = ld(u,htu)

[h, w] = size(u);
dummy = zeros(h,w);

%impose boundary conditions
u = dummies(u);

for i=2:h-1
    for j=2:w-1           
           
        % estimate second-order derivatives
        uxx=u(i,j-1)-2*u(i,j)+u(i,j+1);
        uyy=u(i-1,j)-2*u(i,j)+u(i+1,j); 
        
        % estimate Laplacian
        ddu=uxx+uyy;
        
        x = u(i,j);
           
        % estimate change in the intensity
        del = htu*ddu;
        
        % update intensity
        x = max(0,x+del);
        x = min(255,x);
        dummy(i,j) = x;
    end
end

u = dummy;