clear; clc; close all;

% set the data.  Here I chose to put it in a matrix, but we could use two
% vectors instead...
data = [
   0  4.579
   20 17.535
   40 55.324
   60 149.38
   70 233.7
   80 355.1
   85 433.6
   90 525.76
   95 633.9
   98 707.27
   100 760
   101 787.57 ];

% It will be a bit easier to use T and p rather than
% slicing the data array each time we need it.
T = data(:,1);
p = data(:,2);

% Solve for the coefficients of the polynomial.
% Note that we use a polynomial of order n-1,
% where n is the number of data points.
coefs = polyfit( T, p, length(T)-1 );

% set the points where we want to evaluate the polynomial.  These are the
% points where we want to interpolate to...
x = linspace(0,101);

% evaluate the polynomial at the points given above.
y = polyval( coefs, x );

% plot the original data as well as our polynomial interpolant.  Note that
% we get a lot of "wiggles" in the data.  This is because we are using a
% very high order polynomial.
plot(T,p,'ro',x,y,'k-');