x=linspace(0,0.3,4);
x0=0.15;
y=cos(x);
d=Newtons_Interp(x,y);
dx=x0-x;
y0=d(end);
for i=3:-1:1
  y0=y0*dx(i)+d(i);
end;
printf('y0=%f, cos(%0.2f)=%f\n',y0, x0, cos(x0));

x=linspace(0,0.3,4)';
y=cos(x);
x0=0.15;
A=[ones(size(x)),x];
for i=2:3
  A=[A,A(:,end).*x];
end;
a=A\y;
y0=d(end);
for i=3:-1:1
   y0=y0*x0+a(i);
end;
printf('y0=%f, cos(%0.2f)=%f\n',y0, x0, cos(x0));

x=linspace(0,0.3,4)';
y=cos(x);
x0=0.15;
d=polyfit(x,y,3);
y0=polyval(d,x0);
printf('y0=%f, cos(%0.2f)=%f\n',y0, x0, cos(x0));

xx=linspace(0,0.3);
yp=polyval(d,xx);
yy=cos(xx);
plot(xx,yp,xx,yy,x,y,'o');
figure
plot(xx,yp-yy);

x=linspace(-2,2,5);
y=1./(1+x.^2);
d=polyfit(x,y,4);
xx=linspace(-2,2);
yp=polyval(d,xx);
yy=1./(1+xx.^2);
figure;
plot(xx,yp,xx,yy,x,y,'o');



x=[0;1;2;3;4]; y=[1.5;2.5;3.0;3.5;3.7];
A=[x,ones(size(x))]; 
a=(A'*A)\(A'*y);
printf('a=(%f,%f)\n',a);
a=A\y;
printf('a=(%f,%f)\n',a);
d=polyfit(x,y,1);
printf('d=(%f,%f)\n',d); 
