function [ft] =runimpdiff2d(A,f,ntime,dt)
%
% step in time using implicit solver
% A is the matrix produced by diff2d
% f is the initial image at time t = 0
%
s = size(f);
npix = s(1)*s(2);
ff = reshape(f,npix,1);
ft = zeros(npix,ntime);
ft(:,1) = ff;
K = sparse(eye(npix)) + dt* A;
for k = 2:ntime
    ft(:,k) = K\ft(:,k-1);
end