Yours is a linear problem, however the easiest way to estimate the parameters is likely an unconstrained nonlinear solver, such as fminsearch:
x = ...;
y = ...;
objfcn = @(b,x) b(1).*x.^3 + b(2).*x.^2 + (2*b(1) + 3*b(2)).*x;
[B,resnorm] = fminsearch(@(b) norm(y - objfcn(b,x)), [1;1]);
xv = linspace(min(x), max(x));
figure
plot(x, y, 'pb')
hold on
plot(xv, objfcn(B,xv), '-r')
hold off
grid
