COMMENT Try holistic finite differences on Burgers equation.  
Consider PDEs of the form u_t = u_xx -eps (perturbations) .

Discretise the problem using mixed patch BC.  

Parameter alfa=a/b in my analytic work, treated as small.

(c) Tony Roberts, Mar 2005;

% improve printing
on div; off allfac; on revpri; factor gam,h,eps;

% 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^3=>0, alfa^2=>0 }; 
it:=0$
repeat begin write it:=it+1;
    deq:=df(v,t)-df(v,x,2)
        +eps*(0*v*df(v,x)+c*df(v,x)+b*df(v,x,3)+a*df(v,x,4));
    bcp:=-sub(xi=+r,df(v,x))
        +(1/h)*(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)))))
        +alfa*(-sub(xi=+r,v) +u(j)
            +gam*(+dmu+r/2*ddu)*r
            +gam^2*del2(+dmu+ddu*r/4)*r*(r^2-1)/6
            +gam^3*del2(del2(+dmu+ddu*r/6))*r*(r^2-1)*(r^2-4)/120 
            );
    bcm:=-sub(xi=-r,df(v,x))
        +(1/h)*(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)))))
        -alfa*(-sub(xi=-r,v) +u(j)
            +gam*(-dmu+r/2*ddu)*r
            +gam^2*del2(-dmu+ddu*r/4)*r*(r^2-1)/6
            +gam^3*del2(del2(-dmu+ddu*r/6))*r*(r^2-1)*(r^2-4)/120 
            );
    g:=g+(gd:=(bcp-bcm)/(h*2*r)-intx(deq,xi)/(2*r));
    v:=v+(vd:=h^2*solv(deq+gd,xi)+h*xi*(bcm+bcp)/2);
    showtime;
end until (it>99)or(deq=0)and(bcp=0)and(bcm=0);
ampeq:=sub(xi=0,v)-u(j); % checks amplitude

% 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);

if 1 then begin % 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;

end;