Neumann BC implementation on 2d rectangular grid

Hi,

I have been trying to figure out a solution for this for a long while
and will appriciate any help.

I am trying to implement Finite Differences scheme with Neumann BC to
solve the following PDE:
a*(du/dx) + b*(du/dy) = c
where a,b,c are known.

To do this I have constructed 2d grid, rectangular domain,
for example in matlab:
A=zeros(10,10); A(2:8,2:8)=1;

for points _inside_ the rectangle i.e. A(3:7,3:7) i use central
differences scheme e.g. du/dx = u(x+1,y)-u(x-1,y)/2

the problem occurs with points on the boundary i.e.
A(2:3,:),A(7:8,:),A(:,2:3),A(:,7:8) since i do not know the values on
the boundary, what i did was:
if u(x-1,y) or u(x+1,y) do not exist duplicate in x-direction,
u(x,y-1) or u(x,y+1) do not exist, duplicate in y-direction, in case of
corners i decided to duplicate in the diagonal direction.
And then I can use Neumann BC since the derivative in the normal
direction is 0.
In addition, for the true boundary i use central diff and for the new
duplicated boundary i use forward-backward diff.

This scheme seem to produce errors, what do i do wrong?
can i solve it somehow differently?

Also, there are the points A(1,2),A(2,1),A(8,1),A(1,8),(9,8),(8,9)etc.
that do not have values, i though to give them the average value of
their neigh. but what should be their derivative then?

Thank you very much for your help and time, I'd appriciate any help...
K.

.