function [a,b,c,d,y2,c_flag]=cubsum1(n,x,y,p)% [a,b,c,d,y2,c_flag]=cubsum1(n,x,y,p)% Smoothing cubic spline,% x : x-coordinates monotonic% y : y-coordinates% p : weigthing parameters for each point% a,b,c,d parameters of cubic spline % y2 : second derivatives% dimensions of a,b,c,d,y2, x and y is n% c_flag = 0  normal excecution%        = 1  n<5 not allowed%        = 2  Error in solving system with pentas c_flag = 0; if n<5,   c_flag = 1;   return; end n1=n-1; n2=n-2; for k=1:n1;   d(k)=1.0/(x(k+1)-x(k)); end d(n)=0.0; p1=1.0/p(1); p2=1.0/p(2); h1=d(1); h2=d(2); r1=(y(2)-y(1))*h1; for k=1:n2,   k1=k+1;   k2=k+2;   h3=d(k2);   p3=1.0/p(k2);   s=h1+h2;   a(k)=2.0/h1+2.0/h2+6.0*(h1*h1*p1+s*s*p2+h2*h2*p3);   b(k)=1.0/h2-6.0*h2*(p2*s+p3*(h2+h3));   c(k)=6.0*p3*h2*h3;   r2=(y(k2)-y(k1))*h2;   y2(k)=6.0*(r2-r1);   h1=h2;   h2=h3;   p1=p2;   p2=p3;   r1=r2; end [f_result,c_flag]=pentas(n2,a,b,c,y2); y2=f_result; if c_flag ~=0,    return end for k=n2:-1:1,   y2(k+1)=y2(k); end y2(1)=0.0; y2(n)=0.0; h1=0.0; for k=1:n1,   k1=k+1;   b(k)=d(k);   h2=d(k)*(y2(k1)-y2(k));   d(k)=h2/6.0;   a(k)=y(k)-(h2-h1)/p(k);   c(k)=y2(k)/2.0;   h1=h2; end a(n)=y(n)+h1/p(n); for k=1:n1,   k1=k+1;   h=b(k);   b(k)=(a(k1)-a(k))*h-(y2(k1)+2.0*y2(k))/(6.0*h); endreturn