Now, for a given change in X, Y should change appropriately, however Y is moving in jumps.
For example, instead of moving from 8488 to 8489 as X increases, it waits for X to increase enough, then jumps to 8502.
I know this is some sort of floating point issue (using such small numbers in the polynomial). However I can't put my finger on it, or even see how I can work around it.
XC8 uses reduced 24-bit float format by default which provides insufficient accuracy for your problem. Specify a 32-bit format either for all float variables or specifically for double.
Check if the respective variables are actually created with 32 bit width (e.g. in the memory map or assembly listing).
I was under the assumption that the poor accuracy is due to 24 bit float. But even 32 bit standard IEEE float might be insufficient. One would preferably use 64 bit double for polynomials with large number range. If the x range is limited to a certain positive range (e.g. around 200e3) you can greatly increase the numeric accuracy by rewriting the equations, evolving the polynomial fit around a different origin.