function [M] = relax_G(M,n)
% Relax the interlaced matrix for the G function (using U,V inside M)
% du/dy - dv/dx = G
% Introduce a dummy variable w which represents change in all all
% directions, i.e., solve 
% ((U[i+1,j] + w) - (U[i-1,j] - w))/2h - 
% ((V[i,j+1] - w) - (V[i,j-1] + w))/2h = G[i,j]
%
%
%
s = size(M) ; rows = s(1); cols = s(2);
h = 1/n;

for i= 3:2:rows - 2
    for j = 3:2:cols - 2
        w = .25*(2*h*M(i,j) -(M(i+1,j) - M(i-1,j) - (M(i,j+1) - M(i,j-1))));
        M(i+1,j) = M(i+1,j) + w;
        M(i-1,j) = M(i-1,j) - w;
        M(i,j+1) = M(i,j+1) - w;
        M(i,j-1) = M(i,j-1) + w;
    end
end