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.
%
% Example:
%
% 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, pmanandhar@umassd.edu

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
end
S = B*P;
x = S(:, 1);
y = S(:, 2);

I don't understand the math well enough to pick out exactly what's wrong, but the curve that this code gives for certain sets of points is definitely not right.

Try this and you'll see:
P = [0 2; 4.38 6.38; 5.32 5.32; 9.7 9.7];
x = (0:0.01:9.7);
y = BezierCurve(length(x), P);
plot(x, y, P(:, 1), P(:, 2), 'x-', 'LineWidth', 2); set(gca, 'FontSize', 16)