This is not parameter estimation (‘curve fitting’) as such, but if you want a way to fit all the points exactly, using the interp1 function with the 'spline' method will work:
N = 50; % Number Of Points
x = 1:N;
y = rand(1, N);
xi = linspace(min(x), max(x), 250);
yi = interp1(x, y, xi, 'spline');
figure(1)
plot(x, y, 'gp')
hold on
plot(xi, yi, '-r')
hold off
grid
