improve fit compared to R

Confidential Solution.
Of course it's true, correct, right, proven by facts and sci-tech processes. My results are satisfyying, yup ok. Data analysis: Spline is cubic and as accurate polyfit of order 2, because of the input data, more input will lead to the distinguishing what is more optimal 2 or 3 order of polynomial (spline or quadratic). I have updated, spline and polynomial derivative. It is just left to find derivative of polynomial by only coefficients and not like spline when you find derivative by the full structure(made already) of order in TCE with builtin or workaround functions ... m working on it ... I hope I answered your question.
Here it is the derivative values I found first (for spline and polynome):

dfdx1p = polyder(p2)
dfdx1p =
0.4571 -0.4629

here is the replete, outright answer from me. Please accept.

%Begin:

clc
clear all
close all

%% Measurement(Experiment, Lab) Input Datum:

x_val = [0.858734025700000 0.975248753121875 1.02525125312812 1.13140821289062 1.15927407430000 1.18768630564687 1.21665290240000 1.27628156250000];
y_val = [0.0255926963287944 0.0206749804955220 0.0216695267388959 0.0221177913199036 0.0248476224632830 0.0277388242602486 0.0305936964241147 0.0362101199563917];

%% Real world natural datum to extrapolate to:

x_values_finer = (1+(-0.1))^5:(0.15-(-0.1))/1000:(1+0.15)^5;
%equivalent to:
%x_values_finer = linspace(0.5904900,2.0112400,100);

%% System/function identification by polynomial interpolation/fit:

PP = spline (x_val, y_val)
PP2 = pchip(x_val, y_val)
p2 = polyfit(x_val,y_val,2); % fit a parabola
p22 = mkpp(p2,x_val);
%p22 = polyder(p2);
%p9 = polyfit(x_val,y_val,9);
%% Interpolated system application on real world natural datum(Extrapolation):
YI = spline (x_val, y_val, x_values_finer);
%% In order to compute derivative:
bk = PP.breaks;
cf = PP.coefs;
[m,n] = size(cf);
syms z
v = z.^((n-1:-1:0)');
eq = cf*v;
%% Calculate derivative:
dfdx1 = fnder(PP,1);
dfdx2 = fnder(PP,2);
dfdx1p = polyder(p2);
dfdx2p = polyder(dfdx1p);

%Compute derivative:
der1 = ppval(dfdx1, x_values_finer);
der2 = ppval(dfdx2,x_values_finer);
%Plot derivative:
figure
plot(x_values_finer,der1,'-b',x_values_finer, der2,'-r');
title(' Spline derivative order 1 & 2')

%Compute derivative:
der2p = dfdx2p*ones(1,length(x_values_finer));%= ppval(mkpp(dfdx2p,x_val), x_values_finer);
der1p = x_values_finer*dfdx1p(1) + dfdx1p(2); %ppval(mkpp(dfdx1p,x_val), x_values_finer);

%Plot derivative:
figure
plot(x_values_finer,der1p,'-b',x_values_finer, der2p,'-r');
title('Quadratic poly derivative order 1 & 2')

YI1 = pchip(x_val,y_val,x_values_finer);

f2 = polyval(p2,x_values_finer);
%f9 = polyval(p9,x_values_finer);

%% Validation:
y_int = ppval (PP, x_val);
y_int2 = ppval(PP2,x_val);
y2 = polyval(p2,x_val);
%y9 = polyval(p9,x_val);

%% Resultants:
figure
plot(x_val,y_val)
title('yval vs xval')

figure
plot(x_values_finer)
title('X Values Finer')

figure
plot(x_val,y_int)
title('Y computed from Spline interp.')

figure
plot(YI)
title('Interpolated with Spline values')

figure
plot(x_val,y_int2)
title('Y computed from Piecewise Cubic Hermite Interpolating Polynomial interp.')

figure
plot(YI1)
title('Interpolated with Piecewise Cubic Hermite Interpolating Polynomial values')

figure
plot(x_val,y_val,'o',x_values_finer,f2,'-')
title('Pol1')
%figure
%plot(x_val,y_val,'o',x_values_finer,f9,'*')
%title('Pol2')
figure
plot(x_val,y2,'.')
title('Pol3')
%figure
%plot(x_val,y9,'+')
%title('Pol4')

%Developed from needing help code by
%https://independent.academia.edu/PMazniker
%+380990535261
%https://join.skype.com/invite/oXnJhbgys7oW
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w

%The End.