Hm. I don't know what parula() was originally interpolated in, but I gave it a shot in LAB just for fun. At least it's smooth enough that cubic interpolation behaves fairly well with few points. Of course, I didn't stop to think that lab2rgb() and parula() came out in the same version (R2014b). Still, there are other options.
% a rough trajectory in LAB
x0 = [0 0.07451 0.1294 0.1843 0.2431 0.3647 0.4902 0.5569 0.6039 0.6745 0.7373 0.7843 0.8314 0.8784 0.9412 1];
CT0 = [27 47 -66; 37 53 -80; 43 49 -82; 48 41 -80; 53 26 -73; 63 -9 -44; 69 -38 -12; ...
72 -47 9; 74 -52 26; 75 -41 52; 76 -19 69; 77 0 72; 80 12 68; 84 4 76; 90 -12 84; 96 -22 92];
% generate a new uniform table of desired length
N = 256; % length of new color table
xf = linspace(0,1,N);
CT1 = interp1(x0,CT0,xf,'pchip');
% IPT, R2014b+ (which has parula() anyway)
%CT1 = lab2rgb(CT1);
% IPT, R2010a+ (for WP adaptation option)
cform = makecform('lab2srgb','adaptedwhitepoint',whitepoint('D65'));
CT1 = applycform(CT1,cform);
% MIMT, legacy-friendly (see File Exchange)
%CT1 = ctflop(lch2rgb(lab2lch(ctflop(CT1)),'lab'));
% create test sweeps for comparison
% parula() on top, recreated on bottom
CT = parula(N);
testsweeps = cat(3,CT,CT1);
testsweeps = permute(testsweeps,[3 1 2]);
testsweeps = imresize(testsweeps,[128 256],'nearest');
% compare the results
imshow(testsweeps)
So the difference is noticeable if you have them next to each other, but it's not terrible for a short 16-point reference curve. I might have been able to get away with less in HSYn, but I was trying to stick with base/IPT tools.
EDIT: I went ahead and did it in HSYn. I probably could have placed some of the points better, but it's close enough for arm's length with only 11 points. It doesn't really save any time, but I guess that wasn't the goal. There aren't any base/IPT tools to do this though, you'd need MIMT to run this example.
x0 = [0 0.1647 0.2588 0.3412 0.4471 0.5412 0.6353 0.7608 0.8353 0.9451 1];
CT0 = [252 1.806 0.236; 236 1.488 0.4017; 219 1.376 0.4581; 204 1.217 0.5009; ...
187 1.422 0.5166; 170 1.213 0.5772; 128 1.023 0.6381; 57 1.901 0.7024; ...
41 2.418 0.756; 58 4.209 0.8334; 63 6.642 0.8788];
% generate a new uniform table of desired length
N = 256;
xf = linspace(0,1,N);
CT1 = interp1(x0,CT0,xf,'pchip');
% HSY conversions require MIMT
CT1 = ctflop(hsy2rgb(ctflop(CT1),'native'));
% create test sweeps for comparison
% parula() on top, recreated on bottom
CT = parula(N);
testsweeps = cat(3,CT,CT1);
testsweeps = permute(testsweeps,[3 1 2]);
testsweeps = imresize(testsweeps,[128 256],'nearest');
% compare the results
imshow(testsweeps)
There are other attempts that would work prior to R2014b, but both of the above are closer to the version of parula() that would be around in R2020x.
For comparison, R2017b parula() vs R2014b parula() vs paruly()
