BDF:=proc(k,q)
local p,i,t;
  p:=interp([seq(i*dt,i=0..k)],[seq(y[i],i=0..k)],t);
  simplify(subs(t=q*dt,diff(p,t)))=f(y[q]);
end: