gistfile1.matlab

with(plots);
with(CurveFitting);
spl1 := PolynomialInterpolation([0, 1, 2, 3, 4, 5, 6, 7], [0, 1, 2, 2, 1, 0, 1, 2], t);
spl2 := PolynomialInterpolation([0, 1, 2, 3, 4, 5, 6, 7], [1, 0, 1, 4, 5, 4, 3, 4], t);
plot([spl1, spl2, t=0..7]);
xx := [0, 3, 3, 0]: yy := [4, 9, -5, 0];
n := nops(xx) - 1;
Points := plot([seq([xx[j], yy[j]], j=1..n+1)], style=point);
bez := x->evalm(sum(binomial(n, i)*x^i*(1-x)^(n-i)*[xx[i+1], yy[i+1]], i=0..n));
Pbez := plot([bez(x)[1], bez(x)[2], x=0..1]):
plots[display]([Points, Pbez], axes=framed);