there is a function private/Bspline2pp.m that just does the conversion to pp
It evalutes the spline in k points inside each subitervals then invert the Vandermond matrix to compute a row of pp.coefs. The last stage can be done with polyfit for coding convenience.
You can also compute from Taylor expansion of the polynomial wrt th the left knot as Torsen has suggested.
It requires to compute 0 to (k-1) order derivatives
Code to avoid using fn2fm (why?!) if this idea
WARNING this only works for non-duplicated knot vectors
x = [3.0,4.5,6.0,7.5,9.0,12.0,15.0];
y = [0 0.0343653 0.0694232 0.105143 0.141178 0.246013 0.630537];
f_pp = fn2fm(f_bm, 'pp');
breaks = f_bm.knots(k:end-k+1);
pieces = length(breaks)-1;
coefs(:,k-j) = fnval(Sd, xi)./factorial(j);
pp = struct('form', 'pp', ...
f_pp.coefs
ans =
-0.000010839458734 0.000095398958747 -0.000104662728187 0.022889129608328 -0.000000000000000
0.000030128251860 -0.000067192922266 0.000053996227019 0.023842354392321 0.087249675574952
0.000176086472898 0.000158768966683 0.000311553851941 0.024130562157454 0.132073372169513
pp.coefs
ans =
-0.000010839458734 0.000095398958747 -0.000104662728187 0.022889129608328 -0.000000000000000
0.000030128251860 -0.000067192922266 0.000053996227019 0.023842354392321 0.087249675574952
0.000176086472898 0.000158768966683 0.000311553851941 0.024130562157453 0.132073372169513