clear clc syms x0 x1 x2 x3 h x syms y0 y1 y2 y3 h0=h;%x1-x0; h1=h;%x2-x1; h2=h;%x3-x2; d0=(y1-y0)/(x1-x0); d1=(y2-y1)/(x2-x1); d2=(y3-y2)/(x3-x2); d0=(y1-y0)/h; d1=(y2-y1)/h; d2=(y3-y2)/h; A=[1 0 0 0; h 4*h h 0; 0 h 4*h h; 0 0 0 1 ] b=[0; 6*(d1-d0); 6*(d2-d1); 0]; m=A\b; m0=m(1); m1=m(2); m2=m(3); m3=m(4); s10=y1; s11=d1-h/6*(2*m1+m2); s12=m1/2; s13=(m2-m1)/6/h; s=s10+s11*(x-x1)+s12*(x-x1)^2+s13*(x-x1)^3; ss=s; s=ss; s=subs(s,y1,0); s=subs(s,y2,0); s=subs(s,y3,0); cc=collect(s,y0); simplify(cc); pretty(cc) s=ss; s=subs(s,y0,0); s=subs(s,y2,0); s=subs(s,y3,0); cc=collect(s,y1); simplify(cc); pretty(cc) s=ss; s=subs(s,y0,0); s=subs(s,y1,0); s=subs(s,y3,0); cc=collect(s,y2); simplify(cc); pretty(cc) s=ss; s=subs(s,y0,0); s=subs(s,y1,0); s=subs(s,y2,0); cc=collect(s,y3); simplify(cc); pretty(cc)