General Equation:

a(x,y)(((∂^2)u)/(∂(x^2)))+b(x,y)(((∂^2)u)/(∂(y^2)))+c(x,y),((∂u)/(∂x))+d(x,y)((∂u)/(∂y))+f(x,y)u(x,y)=g(x,y)
(Subscript[h, l])(y)u((Subscript[x, min]),y)+(Subscript[p, l])(y)(Subscript[u, x])((Subscript[x, min]),y)=(Subscript[q, l])(y) 
(Subscript[h, r])(y)u((Subscript[x, max]),y)+(Subscript[p, r])(y)(Subscript[u, x])((Subscript[x, max]),y)=(Subscript[q, r])(y) 
(Subscript[h, b])(x)u(x,(Subscript[y, min]))+(Subscript[p, b])(x)(Subscript[u, y])(x,(Subscript[y, min]))=(Subscript[q, b])(x) 
(Subscript[h, tp])(x)u(x,(Subscript[y, max]))+(Subscript[p, tp])(x)(Subscript[u, y])(x,(Subscript[y, max]))=(Subscript[q, tp])(x)

Our Example:

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))-(1+x,,y),((∂u)/(∂x))+(x-y),((∂u)/(∂y))=0
u(-1,y)=sin(π,y)   (Subscript[u, x])(1,y)=sin(π,y)   u(x,-1)=1-cos(2π,x)    (Subscript[u, y])(x,1)=cos(π,x)

Define Functions

(*Make the domain into a grid*)
xmin=-1;
xmax=1;
ymin=-1;
ymax=1;
xdivisions=100;
ydivisions=200;
dx=((xmax-xmin)/xdivisions);
dy=((ymax-ymin)/ydivisions);
xgrid=Range[xmin,xmax,dx];
ygrid=Range[ymin,ymax,dy];
grid=Outer[{#1,#2}&,xgrid,ygrid];
(*Make an array of values for the solution. Each entry in the array corresponds to the value of the unknown function at a point in space*)
U=Array[u,{xdivisions+1,ydivisions+1},{{xmin,xmax},{ymin,ymax}}];

(*specify functions in above equation*)
a[x_,y_]:=1;
b[x_,y_]:=1;
c[x_,y_]:=-(1+x*y);
d[x_,y_]:=x-y;
f[x_,y_]:=0;
g[x_,y_]:=0;
hl[y_]:=1;
hr[y_]:=0;
hb[x_]:=1;
ht[x_]:=0;
pl[y_]:=0;
pr[y_]:=1;
pb[x_]:=0;
pt[x_]:=1;
ql[y_]:=Sin[π*y];
qr[y_]:=Sin[π*y];
qb[x_]:=1-Cos[2*π*x];
qt[x_]:=Cos[π*x];
(*Create tables of the function values, evaluated at grid points*)
aGrid=Apply[a,grid,{2}];
bGrid=Apply[b,grid,{2}];
cGrid=Apply[c,grid,{2}];
dGrid=Apply[d,grid,{2}];
fGrid=Apply[f,grid,{2}];
gGrid=Apply[g,grid,{2}];
(*Find finite difference appoximations for different derivatives of u*)
dudx=NDSolve`FiniteDifferenceDerivative[{1,0},{xgrid,ygrid},U];
dudx2=NDSolve`FiniteDifferenceDerivative[{2,0},{xgrid,ygrid},U];
dudy=NDSolve`FiniteDifferenceDerivative[{0,1},{xgrid,ygrid},U];
dudy2=NDSolve`FiniteDifferenceDerivative[{0,2},{xgrid,ygrid},U];

Boundary Conditions

(*Consider the left-side bc u(xmax,y)-h1[y]=0. Evaluate the LHS of this equation at each point along the left-side boundary*)
leftbc=Table[hl[ymin+k*dy]*U[[1,k+1]]+pl[ymin+k*dy]*dudx[[1,k+1]]-ql[ymin+k*dy],{k,0,Length[ygrid]-1}];

(*Do the same for the right-side bc*)
rightbc=Table[hr[ymin+k*dy]*U[[-1,k+1]]+pr[ymin+k*dy]*dudx[[-1,k+1]]-qr[ymin+k*dy],{k,0,Length[ygrid]-1}];

(*And for the bottom side boundary condition*)
bottombc=Table[hb[xmin+k*dx]*U[[k+1,1]]+pb[xmin+k*dx]*dudy[[k+1,1]]-qb[xmin+k*dx],{k,0,Length[xgrid]-1}];

(*And for the top side bc*)
topbc=Table[ht[xmin+k*dx]*U[[k+1,-1]]+pt[xmin+k*dx]*dudy[[k+1,-1]]-qt[xmin+k*dx],{k,0,Length[xgrid]-1}];

(*takes the left boundary conditions, and solves them to yield the values of u along the left side *)uleft=NSolve[Map[#==0&,leftbc],Table[U[[1,k]],{k,1,Length[ygrid]}]];

(*takes the right boundary conditions, and solves them to yield the values of u along the right side *)uright=NSolve[Map[#==0&,rightbc],Table[U[[-1,k]],{k,1,Length[ygrid]}]];

(*takes the bottom boundary conditions, and solves them to yield the values of u along the bottom side *)ubottom=NSolve[Map[#==0&,bottombc],Table[U[[k,1]],{k,1,Length[xgrid]}]];

(*takes the top boundary conditions, and solves them to yield the values of u along the top side *)utop=NSolve[Map[#==0&,topbc],Table[U[[k,-1]],{k,1,Length[xgrid]}]];

(*To amke sure boundary conditions are consitent, we only want one u value per point*)
uleft1=uleft[[1,1;;Length[ygrid]]];
utop1=utop[[1,2;;Length[xgrid]]];
uright1=uright[[1,1;;Length[ygrid]-1]];
ubottom1=ubottom[[1,2;;Length[xgrid]-1]];
(*This is a list of all the boundary values of u*)
boundary=Join[uleft1,uright1,utop1,ubottom1];

Solving

(*We now create a set of equations using the PDE. We make a table, with each entry corresponding to an interior grid point. Each entry in the table becomes an equation, from the discretized PDE. We use our knowledge of the boundary conditions to eliminate the the values of u on the boundary*)equations=Map[(0==#)&,Flatten[Table[aGrid[[i,j]]*dudx2[[i,j]]+bGrid[[i,j]]*dudy2[[i,j]]+cGrid[[i,j]]*dudx[[i,j]]+dGrid[[i,j]]*dudy[[i,j]]+fGrid[[i,j]]*U[[i,j]]+gGrid[[i,j]],{i,2,Length[xgrid]-1},{j,2,Length[ygrid]-1}],1]/.boundary];

(*We solve these equations, for the interior values of u*)
intSol=NSolve[equations,Flatten[U[[2;;Length[xgrid]-1,2;;Length[ygrid]-1]]]][[1]];

(*We substitute the interior values of u into the equations that determine the boundary values of u*)
boundarySol=(boundary/.intSol);

(*solutionArray is an array of the values of u. I'm not sure why I had to do the last replacement rule twice!*)
solutionArray=((U/.intSol)/.boundarySol)/.boundarySol;

(*Use this to make a table of 3d coordinates*)
dataPoints=Table[{xmin+i*dx,ymin+j*dy,solutionArray[[i+1,j+1]]},{i,0,Length[xgrid]-1},{j,0,Length[ygrid]-1}];

(*Plot the solution*)
ListPlot3D[Flatten[dataPoints,1],AxesLabel->{Style[x,Medium,Blue],Style[y,Medium,Blue],Style[u,Medium,Blue]},PlotRange->All]