source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
# experimental points that I know have polynomial behavior
x = c(0.8,  2,  3.5,4.9,  6.2,  7.8,  9.3, 11.1,12.5,14.2)
y = c(54.104,59,104,218.3,416.3,821.6,1404,2421,3493,5177)
BF=3.5; HF=2.5
Plane(0,15, 50,5200)
Point(x,y,"brown"); polyli(x,y,"brown")
# I draw the graph of the slope (I take the value at the center
# of each interval)
length(x)
# 10   (they are 10 points)
n = 1:(length(x)-1)
dy = (y[n+1]-y[n])/(x[n+1]-x[n]); dx = (x[n+1]+x[n])/2
# I calculate in which range the slope varies
min(dy); max(dy)
#  4.08   990.5882
# I memorize the previous graph, because then I'll call it
fig = recordPlot()
Plane(0,15, 0,1000)
Point(dx,dy,"brown"); polyli(dx,dy,"brown")
  
# I draw the graph of the slope of the slope
n = 1:(length(Px)-2)
ddy = (dy[n+1]-dy[n])/(dx[n+1]-dx[n]); ddx = (dx[n+1]+dx[n])/2
min(ddy); max(ddy)
#  19.2  145.08
Plane(0,15, 0,150)
Point(ddx,ddy,"brown"); polyli(ddx,ddy,"brown")
# With 2 steps to the slope I arrive at a (almost) linear trend; therefore the
# tabulated function was a polynomial of degree 3
# I draw the graph of the slope of the slope of the slope
n = 1:(length(Px)-3)
dddy = (ddy[n+1]-ddy[n])/(ddx[n+1]-ddx[n]); dddx = (ddx[n+1]+ddx[n])/2
min(dddy); max(dddy)
#  11.28327  12.52744
Plane(0,15, 0,14)
Point(dddx,dddy,"brown"); polyli(dddx,dddy,"brown")
  
# From the fact that the last graph is, approximately, the straight line of ordinate 12 I
# deduce that the leading coefficient of the polynomial is, approximately, 12/3/2 = 2
#
# We verify what we have found by searching directly (with the cubic regression) for
# the 3rd degree polynomial that best approximates the initial data
regression3(x,y)
# 2 * x^3 + -2.99 * x^2 + -0.0629 * x + 55.1 
f = function(x) 2 * x^3 + -2.99 * x^2 + -0.0629 * x + 55.1
# Let's try to superimpose the graph of f on the initial points
replayPlot(fig)
graph1(f, 0,15, "blue")
                       
# OK