x=[4281 6184 6696 8386 8610 10124 10665 11941];
y=[33.4445 42.0849 47.2552 51.0583 55.6392 58.3847 61.8866 64.3721];
x=log(x);y=log(y);
fun=inline('c(1)+c(2)*x','c','x');
c=lsqcurvefit(fun,[5,20],x,y),disp('logy=c(1)+c(2)*logx')
c = -1.8819 0.6465
syms logy logx
logy=vpa(c(1),5)+vpa(c(2),5)*logx