It is actually nonlinear, because it is nonlinear in the parameters (the partial derivatives of the function with respect to each parameter are functions of themselves or of other parameters).
Using lsqcurvefit (Optimization Toolbox), the code is:
n = length(y);
yd = y(2:n);
yi = y(1:n-1);
f = @(B,y) (1 + B(1)).*y - (B(1)./B(2)).*y.^2; % Objective function
yp = linspace(min(yi), max(yi)); % Plot vector for ‘yi’
bc = lsqcurvefit(f, rand(2,1)*100, yi, yd, [1; 1]*eps, [Inf; Inf])
figure(1)
plot(yd, yi, 'xb')
hold on
plot(yp, f(bc,yp), '-r')
hold off
grid
The fit produces:
b = 2.989459285195542E-01
c = 3.990401160101856E-02
b/c = 7.491625942488537E+00
