x = [ 1, 2,  3,  4,  5,   6,   7,   8,   9,  10];
y = [18, 5, 17, 38, 40, 106, 188, 234, 344, 484];

b0 = rand(1,4); % initial guess at best fitting parameters b

optimfun = @(b) curvecost(b,x,y); % optimization function in terms of b

b = fminsearch(optimfun, b0); % optimize!

plot(x, y, 'b.', 'markersize', 30)
hold on

xhat = linspace(min(x),max(x),200);
yhat = curvefun(b,xhat);
plot(xhat, yhat, 'r-', 'linewidth', 2)

function out = curvecost(b,x,y)
    y_est = curvefun(b,x);
    out = sum((y_est-y).^2);

function out = curvefun(b,x)
    out = b(1) + b(2).*x + b(3).*(x.^2) + b(4).*(x.^3);