<PRE class=m name="code">% 分段二次Lagrange插值   
% 命令格式：y=lagrange2(x0,y0,x)   
% x0為節點向量，y0為對應的函數值向量，   
% x為插值點向量，返回值y為x處的函數近似值向量。   
function y=lagrange2(x0,y0,x)   
n=length(x0);   
m=length(x);   
for i=1:m   
   u=x(i);   
   if u<=x0(2)   
      y(i)=y0(1)*(u-x0(2))*(u-x0(3))/(x0(1)-x0(2))/(x0(1)-x0(3))+...   
         y0(2)*(u-x0(1))*(u-x0(3))/(x0(2)-x0(1))/(x0(2)-x0(3))+...   
         y0(3)*(u-x0(1))*(u-x0(2))/(x0(3)-x0(1))/(x0(3)-x0(2));   
   elseif x(i)>x0(n-1)   
      y(i)=y0(n-2)*(u-x0(n-1))*(u-x0(n))/(x0(n-2)-x0(n-1))/(x0(n-2)-x0(n))+...   
         y0(n-1)*(u-x0(n-2))*(u-x0(n))/(x0(n-1)-x0(n-2))/(x0(n-1)-x0(n))+...   
         y0(n)*(u-x0(n-1))*(u-x0(n-2))/(x0(n)-x0(n-1))/(x0(n)-x0(n-2));   
   else   
      for k=2:n-1   
         if u>(x0(k)+x0(k-1))/2 & u<=(x0(k)+x0(k+1))/2   
            y(i)=y0(k-1)*(u-x0(k))*(u-x0(k+1))/(x0(k-1)-x0(k))/(x0(k-1)-x0(k+1))+...   
               y0(k)*(u-x0(k-1))*(u-x0(k+1))/(x0(k)-x0(k-1))/(x0(k)-x0(k+1))+...   
               y0(k+1)*(u-x0(k-1))*(u-x0(k))/(x0(k+1)-x0(k-1))/(x0(k+1)-x0(k));   
         end   
      end   
   end   
end   
y;