Lucy-Richardson algorithm

I am trying to develop an algorithm for Lucy-Richardson deconvolution, but i am facing a problem in making algorithm. The equation which i am using has been attached with this question.

close all; 
clear all; 
%image preprocessing
y = rgb2gray((imread('football.jpg')));
y = im2double(y); 

%take disk psf
PSF = fspecial('disk', 8);

%convolve image with psf
%or use imfilter w/ 'conv'
yblur = imfilter(y,PSF); 

%plot original image
figure(); 
subplot(2,1,1); imshow(y); title('actual image'); 
%plot unnoisy blurred image
subplot(2,1,2); imshow(yblur); title('blurred image');

%zero pad the psf to match the size of the blurred image
%noisy images are all the same size, thus do not require unique PSF's
newh = zeros(size(yblur));  
psfsize = size(PSF);
newh(1: psfsize(1),1:psfsize(2))= PSF;
H = newh;

%Lucy-Richardson Deconvolution
I=zeros(size(y));
J=zeros(size(y));
for x=0:5
    J{x}=((yblur)./(imfilter(I{x},H)));
    I{x+1}=(imfilter(J{x},-H).*I{x});
end
figure
imshow(I)