When I integrate to get N(t) as:
n = dsolve( diff(N) == a*N / (1 + a*h*N), 'IgnoreAnalyticConstraints', true )
in the Symbolic Math Toolbox, the solution is:
N(t) = [0; lambertw(0, a*h*exp(C3 + a*t))/(a*h)]
and then when I have it create an implicit function for the solution, produces:
Nt = @(C3,a,h,t)[0.0; lambertw(0,a.*h.*exp(C3+a.*t))./(a.*h)];
adding N(0) = N0 as an initial condition:
n = dsolve( diff(N) == a*N / (1 + a*h*N), N(0) == N0, 'IgnoreAnalyticConstraints', true )
and integrating yields:
N(t) = lambertw(0, N0*a*h*exp(a*t)*exp(N0*a*h))/(a*h)
and as an anonymous function:
Nt0 = @(N0,a,h,t) lambertw(0,N0.*a.*h.*exp(a.*t).*exp(N0.*a.*h))./(a.*h);
The lambertw function exists in MATLAB (although I admit I've not heard of it until now).
You need to combine parameters N0, a, and h into one vector for curve-fitting purposes in nlinfit or lsqcurvefit. (I do not have the Curve Fitting Toolbox, and use the Statistics and Optimization Toolboxes functions.)
This at least solves your problem of having y on both sides of the equation you want to fit.