General Equation:

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

Our Example:

2
∂
u
∂
2
x
+
2
∂
u
∂
2
y
-(1+xy)
∂u
∂x
+(x-y)
∂u
∂y
=0​​u(-1,y)=sin(πy)
u
x
(1,y)=sin(πy)u(x,-1)=1-cos(2πx)
u
y
(x,1)=cos(πx)

Define Functions
​

In[]:=
(*Makethedomainintoagrid*)​​ClearAll["Global`*"];​​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];​​(*Makeanarrayofvaluesforthesolution.Eachentryinthearraycorrespondstothevalueoftheunknownfunctionatapointinspace*)​​U=Array[u,{xdivisions+1,ydivisions+1},{{xmin,xmax},{ymin,ymax}}];
(*specifyfunctionsinaboveequation*)​​a[x_,y_]:=-(x*y+1);​​b[x_,y_]:=(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];​​(*Createtablesofthefunctionvalues,evaluatedatgridpoints*)​​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}];​​(*Findfinitedifferenceappoximationsfordifferentderivativesofu*)​​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];​​dudxdy=NDSolve`FiniteDifferenceDerivative[{0,2},{xgrid,ygrid},U];
Power
:Infinite expression
1
0
encountered.
Power
:Infinite expression
1
0
encountered.
Power
:Infinite expression
1
0
encountered.
General
:Further output of Power::infy will be suppressed during this calculation.
Power
:Infinite expression
1
0
encountered.
Power
:Infinite expression
1
0
encountered.
Power
:Infinite expression
1
0
encountered.
General
:Further output of Power::infy will be suppressed during this calculation.

Boundary Conditions

In[]:=
(*Considertheleft-sidebcu(xmax,y)-h1[y]=0.EvaluatetheLHSofthisequationateachpointalongtheleft-sideboundary*)​​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}];
In[]:=
(*Dothesamefortheright-sidebc*)​​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}];
In[]:=
(*Andforthebottomsideboundarycondition*)​​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}];
In[]:=
(*Andforthetopsidebc*)​​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}];
In[]:=
(*takestheleftboundaryconditions,andsolvesthemtoyieldthevaluesofualongtheleftside*)uleft=NSolve[Map[#0&,leftbc],Table[U[[1,k]],{k,1,Length[ygrid]}]];
In[]:=
(*takestherightboundaryconditions,andsolvesthemtoyieldthevaluesofualongtherightside*)uright=NSolve[Map[#0&,rightbc],Table[U[[-1,k]],{k,1,Length[ygrid]}]];
In[]:=
(*takesthebottomboundaryconditions,andsolvesthemtoyieldthevaluesofualongthebottomside*)ubottom=NSolve[Map[#0&,bottombc],Table[U[[k,1]],{k,1,Length[xgrid]}]];
In[]:=
(*takesthetopboundaryconditions,andsolvesthemtoyieldthevaluesofualongthetopside*)utop=NSolve[Map[#0&,topbc],Table[U[[k,-1]],{k,1,Length[xgrid]}]];
In[]:=
(*Toamkesureboundaryconditionsareconsitent,weonlywantoneuvalueperpoint*)​​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]];​​(*Thisisalistofalltheboundaryvaluesofu*)​​boundary=Join[uleft1,uright1,utop1,ubottom1];

Solving

In[]:=
(*WenowcreateasetofequationsusingthePDE.Wemakeatable,witheachentrycorrespondingtoaninteriorgridpoint.Eachentryinthetablebecomesanequation,fromthediscretizedPDE.Weuseourknowledgeoftheboundaryconditionstoeliminatethethevaluesofuontheboundary*)equations=Map[(0#)&,Flatten[Table[aGrid[[i,j]]*dudx2[[i,j]]+bGrid[[i,j]]*dudy2[[i,j]]+cGrid[[i,j]]*dudxdy[[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];
Infinity
:Indeterminate expression ComplexInfinity+ComplexInfinity-
99-
2500
3
u-
49
50
,
24
25
+
40000
3
u-
49
50
,
97
100
-25000u-
49
50
,
49
50
+
40000
3
u-
49
50
,
99
100
-
2500
3
u-
49
50
,1
2500
-
49
25
25
3
u-
49
50
,
24
25
-
200
3
u-
49
50
,
97
100
+
200
3
u-
49
50
,
99
100
-
25
3
u-
49
50
,1 encountered.
Infinity
:Indeterminate expression ComplexInfinity+ComplexInfinity-
49
625
-
2500
3
u-
24
25
,
47
50
+
40000
3
u-
24
25
,
19
20
-25000u-
24
25
,
24
25
+
40000
3
u-
24
25
,
97
100
-
2500
3
u-
24
25
,
49
50
-
48
25
25
3
u-
24
25
,
47
50
-
200
3
u-
24
25
,
19
20
+
200
3
u-
24
25
,
97
100
-
25
3
u-
24
25
,
49
50
 encountered.
Infinity
:Indeterminate expression ComplexInfinity+ComplexInfinity-
291-
2500
3
u-
47
50
,
23
25
+
40000
3
u-
47
50
,
93
100
-25000u-
47
50
,
47
50
+
40000
3
u-
47
50
,
19
20
-
2500
3
u-
47
50
,
24
25

2500
-
47
25
25
3
u-
47
50
,
23
25
-
200
3
u-
47
50
,
93
100
+
200
3
u-
47
50
,
19
20
-
25
3
u-
47
50
,
24
25
 encountered.