% 1. polynomials are stored as coefficient vectors
% 3*x^3 - 2.23*x^2 - 5.1*x + 9.8
c = [3 -2.23 -5.1 9.8]

c =

    3.0000   -2.2300   -5.1000    9.8000

% with polyval we can evaluate the polynomial
y = polyval(c,0)

y =

    9.8000

z = polyval(c,1)

z =

    5.4700

sum(c)

ans =

    5.4700

% to plot the polynomial we first sample
x = -1:1/1000:+1;
y = polyval(c,x);
plot(x,y)
% to find the roots we use roots:
r = roots(c);
r

r =

  -1.5985          
   1.1709 + 0.8200i
   1.1709 - 0.8200i

format long e
r

r =

    -1.598531030631471e+000                         
     1.170932181982400e+000 +8.200369994858143e-001i
     1.170932181982400e+000 -8.200369994858143e-001i

y = polyval(c,r)

y =

    -5.151434834260726e-014                         
     7.105427357601002e-015 -8.881784197001252e-015i
     7.105427357601002e-015 +8.881784197001252e-015i

% we see from the size of the evaluation
% that the roots are close to working precision
% We can compute the norm of the residual vector:
norm(y)

ans =

    5.396734625749145e-014

% the backward error is measured by seeing how
% close c is to the coefficient vector of the
% polynomial with roots in r
cr = poly(r)

cr =

  Columns 1 through 2

    1.000000000000000e+000   -7.433333333333296e-001

  Columns 3 through 4

   -1.700000000000005e+000    3.266666666666668e+000

format short e
roots(cr)

ans =

 -1.5985e+000              
  1.1709e+000 +8.2004e-001i
  1.1709e+000 -8.2004e-001i

r

r =

 -1.5985e+000              
  1.1709e+000 +8.2004e-001i
  1.1709e+000 -8.2004e-001i

% poly(r) returns the monic polynomial that
% has roots in r, monic = first coefficient = 1
c

c =

  3.0000e+000 -2.2300e+000 -5.1000e+000  9.8000e+000

c1 = c/c(1)

c1 =

  1.0000e+000 -7.4333e-001 -1.7000e+000  3.2667e+000

cr

cr =

  1.0000e+000 -7.4333e-001 -1.7000e+000  3.2667e+000

% the backward error is
norm(c1 - cr)

ans =

  6.1224e-015

% Let us look at a random polynomial
c = rand(1,100);
r = roots(c);
figure
plot(r,'r+')
%
% 2. curve fitting
%
x = -1:1/100:+1;
y = polyval(c,x);
% we took the random polynomial and evaluated
% at 101 sample points in [-1,+1], so now we have
% data stored at (x,y)
c1 = polyfit(x,y,1)

c1 =

  3.2759e+000  1.3215e+000

c2 = polyfit(x,y,2)

c2 =

  5.7057e+000  3.2759e+000 -5.9941e-001

% In c1 we have the coefficients of a line that
% best approximates the give data.
% In c2 we used a quadratic approximation.
% We will plot the data and the two approximations
figure
plot(x,y,'bx')
hold on
plot(x,polyval(c1,x),'r')
plot(x,polyval(c2,x),'g')
diary off