1;

function B = bezier(P, t)
  n = length(P)-1;
  i = 0:n;
  a = bincoeff(n, 0:n) ;
  b = bsxfun(@power, 1-t, (n-i)');
  c = bsxfun(@power,  t,  i');
  d = bsxfun(@times, a', b.*c)';

  B = permute(sum(bsxfun(@times, d, permute(P, [3,2,1])),2),[3,1,2]);
endfunction

% bezier control points
P = [ 0 1 1
      0 0 1 ];
ts = 0:0.1:1;

plot(B(1,:), B(2,:));