(* ::Package:: *)

function q = CubicValues(p,n)
% p is a column 4-vector with the property that
% p(1) = c(x0)
% p(2) = c(x0+h)
% p(3) = c(x0+2h)
% p(4) = c(x0+3h)
%
% where c(x) is some cubic polynomial, x0 is real and h>0 is real.
%
% q is a column n-vector (n>=4) with the property that
% q(i) = c(x0+(i-1)h), i=1:n

% Initialize q...
q = zeros(n,1);
q(1:4) = p;

% The differences...
d1 = p(4); % = c{x0+3h}
p(4) = p(4)-p(3);
p(3) = p(3)-p(2);
p(2) = p(2)-p(1);
d2 = p(4); % = c{x0+2h,x0+3h}
p(4) = p(4)-p(3);
p(3) = p(3)-p(2);
d3 = p(4); % = c{x0+h,x0+2h,x0+3h}
p(4) = p(4)-p(3);
d4 = p(4); % = c{x0,x0+h,x0+2h,x0+3h}

% By updating the differences we generate the cubic evaluations..
for k=5:n
d3 = d3+d4;
d2 = d2+d3;
d1 = d1+d2;
q(k) = d1;
end

% Correct results: 2 points
% If you use multiplication or power incorrectly, you may lose 1~3 points
% If your time complexity is larger than O(deg^2 + n*deg), you will also lose 1~2 points. deg is the degree of the polynomial, which is 3 in this problem.
%
% In addition, you only need p, q, and an array of length deg for the computing.
% A program that uses deg as a parameter is appreciated.
% It would be great if you made this two thing, though there is no penalty for not doing that.