Derive an FD stencil for non-constant-coefficient (eps u')'

kill(all);

stencil(u):=a*u(-h) + b*u(0) + c*u(h);

depends(eps, x);

make_equation(f):=block([rhs],
  stencil_rhs:diff(eps*diff(f, x), x),
  (      
  stencil(lambda([arg_x], subst([x=arg_x], f)))
  =      
  subst([diff(eps, x)=deps, x=0], stencil_rhs)
  ));                
eqns:[       
  make_equation(1),
  make_equation(x),
  make_equation(x**2)
  ];
display(eqns);
coeffs: solve(eqns, [a,b,c]);
display(coeffs);