Poisson solver for background oriented schlieren

Hi, I'm trying to construct the density field from an optical method called Background oriented schlieren: so from two displacement matrix ux and uy i'm creating the rhs of the poisson eq : Nxx+Nyy=K*rhs the result that i'm getting is not what i'm supposed to get, could any one spot the errors in point of view image analysis and poisson solver

 if true
clear all
close all
ux1=load('bos\ux_pixi2_bos');
ux=ux1.ux;
uy1=load('bos\uy_pixi2_bos');
uy=uy1.uy;                                  
 %%camera parameters and fov
coe_n= -2.911e-04%mm^-2
b_fld=0.1233;  %mm
b_brd=0.1983;%mm
M=0.1983;%mm/pixel
ux=M.*ux;%mm
uy=M.*uy;%mm
% create the RHS element
% Preallocation for the variables du and duy;
[width,height]=size(ux);
dux=zeros(width,height);
duy=zeros(width,height);
%%Compute the central differences 
for k=2:width-1
for j=2:height-1
    dux(k,j)=(ux(k+1,j)-ux(k-1,j))/2*b_brd;
    duy(k,j)=(uy(k,j+1)-uy(k,j-1))/2*b_brd;
end
end
rhs=dux+duy;
rhs(isnan(rhs)) = 0 ;
b=(coe_n).*(rhs);
%%Specifying parametersfor poisson solver
[Nx,Ny]=size(rhs);        
%%Initial Conditions
pn=1.1.*ones(Nx,Ny);        %Preallocating pn
i=2:Nx-1;
j=2:Ny-1;
% Poisson equation solution (iterative) method
tol = 1e-20;             % Set tollerance
maxerr = 1;              % initial error
iter = 0;
p=pn;
 while maxerr > tol
        iter = iter + 1;
% Dirichlet BC
            p(1,:)=1.000227;  
            p(:,Ny)=1.000277;
            p(:,1)=1.000227;

        for j = 2:Ny-1
            for i = 2:Nx-1
 p(i,j)=0.25*(double(pn(i+1,j))...
              +double(p(i-1,j))...
              +double(pn(i,j+1))...
              +double(p(i,j-1))-double(b(i,j)*(b_fld^2)));

            end
            %%Neuman boundary condition along the bottom wall for
      p(Nx,j)=1/4*(2*pn(Nx-1,j)+pn(Nx,j+1)+pn(Nx,j-1)-((b_fld^2)*b(Nx,j)));
         end
   maxerr =sqrt(sum(sum((p-pn).^2)));
  disp(['Maximum error is  ',num2str(maxerr),'Iteration no. ',num2str(iter)]);
   % update pn
    pn=double(p);
    end

regards