COMMENT Try holistic finite differences on a generalised Burgers 
equation.  Consider PDEs of the form u_t = u_xx +... .

Discretise the problem using patch IBC a la Kevrekidis.
This is high order consistent for even and odd operators.  

Tony Roberts, Sept 2004;

% improve printing
on div; off allfac; on revpri; factor gam,h,a,b,c;

% make function of xi=(x-x_j)/h
depend xi,x;
let df(xi,x)=>1/h;

% parameterise with evolving u(j)
operator u; depend u,t;
let df(u(~k),t)=>sub(j=k,g);

% intx computes \int_{-r}^{+r} ... dxi
operator intx; linear intx;
let { intx(xi^~p,xi) => (r^(p+1)-(-r)^(p+1))/(p+1)
    , intx(xi,xi) => 0 
    , intx(1,xi) => 2*r };

% solves v_xixi=RHS s.t. v(0)=0 and v_xi(+r)+v_xi(-r)=0
operator solv; linear solv;
let { solv(xi^~p,xi) => (xi^(p+2)/(p+2)-(r^(p+1)+(-r)^(p+1))*xi/2)/(p+1)
    , solv(xi,xi) => (xi^3-3*r^2*xi)/6
    , solv(1,xi) => (xi^2)/2 };

% linear solution in jth interval
v:=u(j);
g:=0;
% bcs expressed in terms of \mu\delta\pm r\delta^2
procedure del2(f); sub(j=j+1,f)-2*f+sub(j=j-1,f);
dmu:=(u(j+1)-u(j-1))/2$
ddu:=del2(u(j))$

% iteration to some order
% eps measures the amplitude of non-diffusive terms
let { gam^4=>0, eps^2=>0 }; 
it:=0$
repeat begin write it:=it+1;
    deq:=df(v,t)-df(v,x,2)+eps*(c*df(v,x)+b*df(v,x,3)+a*df(v,x,4));
    bcp:=-sub(xi=+r,h*df(v,x))
		+gam*(dmu+r*ddu)
		-gam^2*del2(dmu*(1/6-r^2/2)+ddu*r*(1/12-r^2/6)) 
		+gam^3*del2(del2(dmu*(1/30-r^2/8+r^4/24)
			+ddu*r*(1/90-r^2/36+r^4/120))) 
		-gam^4*del2(del2(del2(dmu*(1/140-r^2*7/240+r^4/72-r^6/720)
			+ddu*r*(1/560-r^2*7/1440+r^4/480-r^6/5040))))
        ;
    bcm:=-sub(xi=-r,h*df(v,x))
		+gam*(dmu-r*ddu)
		-gam^2*del2(dmu*(1/6-r^2/2)-ddu*r*(1/12-r^2/6)) 
		+gam^3*del2(del2(dmu*(1/30-r^2/8+r^4/24)
			-ddu*r*(1/90-r^2/36+r^4/120))) 
		-gam^4*del2(del2(del2(dmu*(1/140-r^2*7/240+r^4/72-r^6/720)
			-ddu*r*(1/560-r^2*7/1440+r^4/480-r^6/5040))))
		;
    g:=g+(gd:=(bcp-bcm)/(h^2*2*r)-intx(deq,xi)/(2*r));
    v:=v+(vd:=h^2*solv(deq+gd,xi)+xi*(bcm+bcp)/2);
    showtime;
end until (it>50)or(deq=0)and(bcp=0)and(bcm=0);
ampeq:=sub(xi=0,v)-u(j); % checks amplitude

% remove artificial parameter
eps:=1;

% in equivde_r$
% deduce the equivalent differential form, 
% du(n) denotes nth derivative of u at the point j
operator du; 
o:=10;
let h^~p=>0 when p>o;
ge:=(g where u(~k)=>(du(0)+for n:=1:o+3 sum h^n*(k-j)^n*du(n)/factorial(n)))$
gee:=sub(gam=1,ge)$
coeff(gee,h);

% deduce operator form
let mu^2=>1+delta^2/4;
gop:=(g where {u(j)=>1, 
    u(j+~p)=>(1+mu*delta+delta^2/2)^p when p>0,
	u(j-~p)=>(1-mu*delta+delta^2/2)^p when p>0 });
factor mu,delta;
vop:=(v where {u(j)=>1, 
    u(j+~p)=>(1+mu*delta+delta^2/2)^p when p>0,
	u(j-~p)=>(1-mu*delta+delta^2/2)^p when p>0 })$

end;