Here is modified code that works. Note that the function now returns two sets of variables.
function [x y] = BezierCurve(N, P)
% This function constructs a Bezier curve from given control points. P is a
% vector of control points. N is the number of points to calculate.
% P = [0 0; 1 1; 2 5; 5 -1];
% y = BezierCurve(1000, P);
% plot(x, y, P(:, 1), P(:, 2), 'x-', 'LineWidth', 2); set(gca, 'FontSize', 16)
% Prakash Manandhar, firstname.lastname@example.org
Np = size(P, 1);
u = linspace(0, 1, N);
B = zeros(N, Np);
for i = 1:Np
B(:,i) = nchoosek(Np-1,i-1).*u.^(i-1).*(1-u).^(Np-i); %B is the Bernstein polynomial value
S = B*P;
x = S(:, 1);
y = S(:, 2);