Arduino 10bit ADC is not accurate?
As stated by Majenko, the ADC can in principle have all sorts of
errors. However, if you look at your data, you should notice that the
error is essentially linear with the voltage. This means that the
non-linearity error is quite small, and you are mainly seeing the offset
and gain errors. This is fortunate, because these errors are easy to
remove, given good calibration data.

The table you wrote in in your question happens to be a perfectly good
calibration table. By fitting a straight line to it you can get a
calibration function, and then correct the readings simply by adding
this to your code:

voltage = 0.9848815 * voltage + 0.125354;


Or you could instead figure out the correction at the ADC level and
replace

adc = analog * 1.228 / 1023;


by

adc = (analog + offset) / scale;


where offset and scale are computed from the calibration.



Edit: Here is how I extracted the calibration coefficients from the
data.

First, I saved the data in a text file calibration.tsv and cleaned it up
to remove the units and add a comment marker (#) in front of the column
labels:

#Input | Arduino | Difference
 1      0.89     -110e-3
 2      1.90     -100e-3
...


Then I looked at the data with gnuplot. Since I want to know how to
correct the readings, I plot the amount I have to add (minus column 3)
v.s. the Arduino reading (column 2):

plot 'calibration.tsv' using 2:(-$3)


At this point it appears it's linear with what looks like mostly
quantization noise. I model this as a line of equation a*x+b, and
compare the line with the data:

a = -1e-2
b = 0.12
plot 'calibration.tsv' using 2:(-$3), a*x+b


The line is quite off the data, but it's a reasonable first
approximation. Now I ask gnuplot to perform a least squares fit, then
take another look:

fit a*x+b 'calibration.tsv' using 2:(-$3) via a, b
replot


Now the fit looks quite good. When it completed the fit, gnuplot
displayed:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = -0.0151185       +/- 0.0002953    (1.953%)
b               = 0.125354         +/- 0.002173     (1.734%)


so now I know I can correct the data with:

voltage += -0.0151185 * volatage + 0.125354;


which is equivalent to what I wrote above.