t = [2,4,6,12,16,30,60,110,180,280,400,515];

ans =

7.8251   10.4813   12.3635   16.5756   19.0916   25.3864   32.8043   36.6479   35.9171   29.3685   16.7998    0.5427


function u = cubic( t )
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
x = [0.52,0.52,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468,494,507,520,520];
y = [5.28794,5.28794,9.4,13.84,20.2,24.9,28.44,31.1,35,36.5,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2,0.2];
n = length(x);
m = n-1;
q = x(1,(2:n-1));
p = y(1,(2:n-1));
a = zeros(1,n-1);
a(1) = 1.86548;
a(n-1) = -0.046115;
for i = 2:1:n-2
    c = q(i) - q(i-1);     %f(xi,xi+1)
    d = p(i) - p(i-1);
    a(i) = d/c;
end
b = zeros(1,n-2);
for i = 2:1:n-1
    c = x(i+1) - x(i-1);       %f(xi,xi+1,xi+2)
    d = a(i) - a(i-1);
    b(i-1) = d/c;
end
e = zeros(1,n-2);
f = e;
for i = 2:1:n-3
    c = q(i) - q(i-1);
    d = q(i+1) - q(i-1);
    e(i) = c/d;
    f(i) = 1 - e(i);
end
v = 2*ones(1,n-2);
g = diag(v,0);
g(1,2) = 1;
g(n-2,n-3) = 1;
for i = 2:1:n-3
    g(i,i-1) = e(i);
    g(i,i+1) = f(i);
end
b  = 6*b;
b = b';
h = g\b
s = length(t);
u = zeros(1,s);
for j = 1:1:s
    if t(j)<= q(1)||t(j)>=q(n-2)
        disp('error');
        break;
    end
    for i = 2:1:n-2
        if t(j)>=q(i-1)&&t(j)<=q(i)
            u(j) = p(i-1)+(p(i)-p(i-1))/(q(i)-q(i-1))*(t(j)-q(i-1));
            u(j) = u(j)-(1/6*h(i)+1/3*h(i-1))*(q(i)-q(i-1))*(t(j)-q(i-1));
            u(j) = u(j)+1/2*h(i-1)*(t(j)-q(i-1))^2;
            u(j) = u(j)+1/6*(h(i)-h(i-1))/(q(i)-q(i-1))*(t(j)-q(i-1))^3;
        end
    end
end
end