kill(all);

g(x):=a*x^2+b*x+c;

linsolve(
    [
        g(f-df)=g1,
        g(f   )=g2,
        g(f+df)=g3
    ],[a,b,c]
),globalsolve=true;

ev(g(x))$
diff(%,x)$
ratsimp(%)$
linsolve(%,x)$
fullratsimp(ev(x,%))$
expand(num(%))/denom(%);
xstat:%$

fullratsimp(g(xstat))$
expand(num(%))/denom(%);
gstat:%$