robin boundary condition :  'bc of the third kind '
peter.vlasschaert@gmail.com            ,01/11/2017

kill(all)$

first order derivative:
three point backward differences formula  ( =numerical approximation)

y(x):=a_[0]*x^0+a_[1]*x^1+a_[2]*x^2;

ref : point x=x0 :'point at the boundary'
         gridspacing = h
see also 'p16' below

inside domain : the two other points (2e,3e)
1e)  ref : x0             =>  y0
2e )          x0-h         =>  y1
3e )          x0-2*h     =>  y2

p1:y0=y(x0);
p2:y1=y(x0-h);
p3:y2=y(x0-2*h);

use : menu => Equations => Solve

p4:ev(solve([p1,p2,p3],[a_[0],a_[1],a_[2]]),x0=0)[1];

solution :first order derivative:  'coef : a_[1] =a   '
                                                                                       1

p5:p4[2];

symbolic: ∂u/∂x = a1, u=u(x,y)    "x direction"
y0 = u(i,j)         also 'ref'
y1 = u(i-1,j)
y2 = u(i-2,j)
rem :point i,j grid point from boundary are the inside of the domain.

p6:'diff(u,x) = rhs(p5);

symbolic: ∂u/∂y = a1, u=u(x,y)     "y direction"
z0 = u(i,j)
z1 = u(i,j-1)
z2 = u(i,j-2)
rem :point i,j grid point from boundary are the inside of the domain.

p7:'diff(u,y) = ev(rhs(p5),y0=z0,y1=z1,y2=z2);

robin bc : x direction   'p8'
                    y direction   'p9'
bc = (b)oundary (c)ondition

p8a:y0 = u(i,j);
p8b:y1 = u(i-1,j);
p8c:y2 = u(i-2,j);
p9a:z0 = u(i,j);
p9b:z1 = u(i,j-1);
p9c:z2 = u(i,j-2);

symbolic representation        : ' robin bc '

p8 :lhs(p6) +c_[1]*'u = c_[2];
p9 :lhs(p7) +c_[3]*'u = c_[4];

numerical representation    : ' robin bc '

p10:ev(rhs(p6),p8a,p8b,p8c)+c_[1]*rhs(p8a) = c_[2];
p11:ev(rhs(p7),p9a,p9b,p9c)+c_[3]*rhs(p9a) = c_[4];

robin bc ' x direction '

p10a:radcan(lhs(p10) - rhs(p10) =0);
p10b:part(p10a,1,1)=rhs(p10a);

robin bc ' y direction '

p11a:radcan(lhs(p11) - rhs(p11) =0);
p11b:part(p11a,1,1)=rhs(p11a);

corner node : crossing from ' robin bc ,x direction and robin bc ,y direction '

p12:p10b+p11b$
p13:ratsimp(p12);

we want to solve : u(i,j) = ?   

p14:solve(p13,u(i,j))[1];

u(i,j)  for other respresentation : (index notation)

p14a:uu(i,j):= uu[i,j];
p14b:ev(p14,u=uu)$
p14c:ev(p14b,uu=u);

when corner coordinate (n,m)           =>       u(n,m)
=>                   "  we treat the corner node as an internal node "
rem : internal node = internal gridpoint ( inside domain            = internal )
rem : ghost    node = ghost     gridpoint ( outside the domain  = external)

p15:ev(p14,i=n,j=m);

u(n,m)  for other respresentation : (index notation)  => 'corner node'

p15a:ev(p15,u=uu)$
p15b:ev(p15a,uu=u);

rem: y(x) ,general ( by use yy(x) )
( grid spacing =constant = h  'value')

p16:yy(x):=sum(a_[k]*x^k, k, 0,3);

p17:yya=yy(x0-h);
p18:yyb=yy(x0);
p19:yyc=yy(x0+2*h);
p20:yyd=yy(x0+3*h);