function J = bh_lst(p0,tdata,pdata)
% Least Squares fit to Beverton_Holt Growth
N = length(tdata);
p = p0(1);
pop = [p];
err = [pdata(1) - p];
for i = 2:N  % Beverton_Holt iteration
    p = p0(2)*p/(1 + p/p0(3));
    pop = [pop,p];
    err = [err, pdata(i) - p];
end
J = err*err'; % Sum of square errors
end

