clear all;
close all;
t=0:0.1:3

X0=3;X1=-3;X2=-3;X3=3;
Y0=3;Y1=3;Y2=-3;Y3=-3;
M0=(0-1)*(0-2)*(0-3);
M1=(1-0)*(1-2)*(1-3);
M2=(2-0)*(2-1)*(2-3);
M3=(3-0)*(3-1)*(3-2);
myx=X0*(t-1).*(t-2).*(t-3)/M0+X1*(t-0).*(t-2).*(t-3)/M1+X2*(t-0).*(t-1).*(t-3)/M2+X3*(t-0).*(t-1).*(t-2)/M3;
myy=Y0*(t-1).*(t-2).*(t-3)/M0+Y1*(t-0).*(t-2).*(t-3)/M1+Y2*(t-0).*(t-1).*(t-3)/M2+Y3*(t-0).*(t-1).*(t-2)/M3;
plot(myx, myy);
hold on;
plot([3, -3, -3, 3],[3, 3, -3, -3],'*')
plot([X0, X1, X2, X3],[Y0, Y1, Y2, Y3],'*')
figure

t=0:0.1:4
X0=3;X1=0;X2=1;X3=0;X4=3;
Y0=2.05;Y1=3;Y2=2;Y3=1;Y4=1.95;
M0=      (0-1)*(0-2)*(0-3)*(0-4);
M1=(1-0)      *(1-2)*(1-3)*(1-4);
M2=(2-0)*(2-1)      *(2-3)*(2-4);
M3=(3-0)*(3-1)*(3-2)      *(3-4);
M4=(4-0)*(4-1)*(4-2)*(4-3)      ;
myx= X0*(t-1).*(t-2).*(t-3).*(t-4)/M0 ...
   + X1*(t-0).*(t-2).*(t-3).*(t-4)/M1 ...
   + X2*(t-0).*(t-1).*(t-3).*(t-4)/M2 ...
   + X3*(t-0).*(t-1).*(t-2).*(t-4)/M3 ...
   + X4*(t-0).*(t-1).*(t-2).*(t-3)/M4;
myy= Y0*(t-1).*(t-2).*(t-3).*(t-4)/M0 ...
   + Y1*(t-0).*(t-2).*(t-3).*(t-4)/M1 ...
   + Y2*(t-0).*(t-1).*(t-3).*(t-4)/M2 ...
   + Y3*(t-0).*(t-1).*(t-2).*(t-4)/M3 ...
   + Y4*(t-0).*(t-1).*(t-2).*(t-3)/M4;
plot(myx, myy);
hold on;
%(3,2)->(0,3)->(1,2)->(0,1)->(3,2)
plot([X0, X1, X2, X3, X4],[Y0, Y1, Y2, Y3, Y4],'*')
plot(myx, myy);