function yi=lsqint1(x0,y0,x1,y1,x2,y2,x3,y3,xi)
% interpolate using least squares quadratic polynomial method
x(1)=x0;	y(1)=y0;
x(2)=x1;	y(2)=y1;
x(3)=x2;	y(3)=y2;
x(4)=x3;	y(4)=y3;
p=polyfit(x,y,2);
yi=polyval(p,xi);