Question: Finite difference Method

ode:=diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)=0

bc1:=df(xmin)=2*sigma__1;

bc2:=d2f(xmin)=0;

bc3:=d3f(xmax)=2;

dfde:=proc(h)(f[k+1]-f[k-1])/2/h;end proc:

dfde(h);

d2fde2:=proc(h)(f[k+1]-2*f[k]+f[k-1])/h^2;end proc:

d2fde2(h);

d3fde3:=proc(h)(f[k+1]-3*f[k]+3*f[k-1]-f[k-2])/h^3;end proc:

d3fde3(h);

d2fde2f:=proc(h)(f[k+2]-2*f[k+1]+f[k])/h^2;end proc:

d2fde2f(h);

d2fde2b:=proc(h)(f[k]-2*f[k-1]+f[k-2])/h^2;end proc:

d2fde2b(h);

dfdef:=proc(h)(f[k+1]-f[k])/h;end proc:

dfdef(h);

h:=xmax/(n-1)

stencil:=subs(diff(f(eta),eta$3)=d3fde3(h),f(eta)=f[k],diff(f(eta),eta$2)=d2fde2,ode);

bcEqs:=[subs(k=0,dfdef(h))=rhs(bc1),subs(k=0,d2fde2f(h))=rhs(bc2),
subs(k=n-1,d2fde2b(h))=rhs(bc3)];

eqs:=Vector(n-2):
cnt:=0:

for k from 1 to n-2 do
    cnt:=cnt+1:
    eqs(cnt):=stencil;
end do:

eqs: