% cubic.m fits a cubic polynomial to two data points

D = [[2 1]; [5 9]];                 % list of data points
x = D(:,1);                         % x values are 1st column of D
y = D(:,2);                         % y values are 2nd column of D

A = [[ones(size(x)) x x.^2 x.^3]; [zeros(size(x)) ones(size(x)) 2*x 3*x.^2]];
                                    % coefficient matrix

c = [y; zeros(size(y))];            % right hand side of Ab = c

b = A\c;                            % solve the equation Ab = c for b

xi = min(x):0.1:max(x);             % x values for interpolation
yi = b(1) + b(2)*xi + b(3)*xi.^2 + b(4)*xi.^3; % interpolated y values

plot(x, y, 'o', xi, yi);            % plot data points and interpolated curve
title('Cubic with horizontal tangents');
axis([0 max(x)+1 0 max(y)+1]);      % choose nice axes

b                                   % display the coefficients