kill(all)$

maxima:programming,'automatic generation of finite difference',version 1.0
peter.vlasschaert@gmail.com,16/05/2020.

load(lrats)$

"Δ2f"(f):=subst(x+1,x,f)-2*f+subst(x-1,x,f);

Δ2f(x^2);

"Δ2fxy"(f,a,b):=lratsubst([x=x+a,y=y+b],f) -2*f+lratsubst([x=x-a,y=y-b],f);

Δ2fxy(f(x,y),h,k);

Δ2fxy(f(x^2),0,1);

"Δ2fxyi"(f,x,y,a,b,i):=lratsubst([x=x+a,y=y+b],f[i+1](x,y))-2*f[i](x,y)+lratsubst([x=x-a,y=y-b],f[i-1](x,y));

Δ2fxyi(f,x,y,a,b,i);

Δ2fxyi(g,x,y,h,k,j);

eq[j] :=Δ2fxyi(u,x,y,h,k,j);

diff(u[0](x,y),x,2) =lratsubst([k=0],eq[1])/h^2;

p2:diff(u[0](x,y),y,2) =lratsubst([h=0],eq[1])/k^2;

q[i,h]:=diff(u[i](x,y),x,2) =lratsubst([k=0],eq[i])/h^2;

m:[];

for i:1 thru 2 do(
m:endcons(q[i,h],m)
);

m;

eqq[j] :=Δ2fxyi(u,ii,jj,h,k,j);

eqq[1];

;qq[i,h]:='diff(u[i](ii,jj),x,2) =lratsubst([k=0],eqq[i])/h^2;

mm:[];

for i:0 thru 2 do(
mm:endcons(qq[i,h],mm)
);

mm;

mmm:ev(mm,h=1);