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)