Image Deblurring using Conjugate Gradient Method

I am working on debluring and image using the conjugate gradient method and I can seem to make any progress. I get a big issue that I have is that I end up attempting to divide scalar my matrices based on the formulas that I have. If anyone could have a look and lend a hand that'd be great. Or if anyone is aware of a sample of this online that I could take a look at?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT u0 = noisy image 
% OUTPUT: I(i,j) = denoised image 

% read the noisy data I0 in .jpg format

Io=imread('figure.jpg'); %opens image and assigns variable Io
Io=double(Io); %returns double precision value for u0
I=Io; %sets u to equal u0 so that u will give a new figure after H1 regularization
[M,N]=size(I); %sets size to equal I
imagesc(Io); axis image; axis off; colormap(gray);

%         PARAMETERS  %
%-----------------------------------------------

% coefficient of the H1 norm (needs to be adapted for each image) 
lambda=0.028; 

%  space discretization 
h=1.;

%  number of iterations (depends on the image) 
IterMax=100;

%  needed to regularize at the origin 
eps=0.000001;

%-----------------------------------------------------------
%---------------------------------------------------------------------------

b = lambda*Io;

Step 1:

%-----------------------------------------------------------
 %%initial guess

for i = 1:length(M)
 for i=2:M-1,
      for j=2:N-1, 

        r(i,j) =  - (I(i,j+1)-2*I(i,j)+I(i,j-1)-(I(i+1,j)-2*I(i,j)+I(i-1,j))); 
        ep = ones(253,252)*eps;

        end
   end
   end

Step 2:

         mu = norm(r)*norm(r); %step two
         r = ep + r;
         p = r;

Step 3:

%-----------------------------------------------------------
%for Iter=1:IterMax, %1 to 100 for loop

  %Iter     

for Iter=1:IterMax, %1 to 100 for loop

  Iter

      for i=2:M-1
        for j=2:N-1 

       q(i,j) = lambda*I(i,j) -(I(i,j+1)-2*I(i,j)+I(i,j-1))-(I(i+1,j)-2*I(i,j)+I(i-1,j));
       q(i,j) = q(i,j) + eps;
       p(i,j) = p(i,j) + eps;
        end
      end

     for i=2:M-1
        for j=2:N-1

       alpha =mu./(p(i,j)'*q(i,j)); %(3) alpha = mu/(p*q)        
       I(i,j) = I(i,j) + alpha * p(i,j); %
       r(i,j) = r(i,j) - alpha * q(i,j);
       mu_pre = mu;
       mu = norm(r)*norm(r);
       p(i,j) = r(i,j) + mu/mu_pre*p(i,j);
       I(i,j) = I(i,j)-alpha*p(i,j);    

        end
     end

%end
%end

figure 
imagesc(I); axis image; axis off; colormap(gray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%