function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]=
leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})

Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).

Version 3.beta
Optional parameters are in braces {}.
x = column vector or matrix of independent variables, 1 row per
observation: x = [x0 x1....xm].
y = column vector of observed values, same number of rows as x.
wt = column vector (dim=length(x)) of statistical weights.  These
should be set to be proportional to (sqrt of var(y))^-1; (That is,
the covariance matrix of the data is assumed to be proportional to
diagonal with diagonal equal to (wt.^2)^-1.  The constant of
proportionality will be estimated.); default = ones(length(y),1).
pin = column vec of initial parameters to be adjusted by leasqr.
dp = fractional increment of p for numerical partial derivatives;
default = .001*ones(size(pin))
dp(j) > 0 means central differences on j-th parameter p(j).
dp(j) < 0 means one-sided differences on j-th parameter p(j).
dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
F = name of function in quotes; the function shall be of the form y=f(x,p),
with y, x, p of the form y, x, pin as described above.
dFdp = name of partial derivative function in quotes; default is "dfdp", a
slow but general partial derivatives function; the function shall be
of the form prt=dfdp(x,f,p,dp,F) (see dfdp.m).
stol = scalar tolerance on fractional improvement in scalar sum of
squares = sum((wt.*(y-f))^2); default stol = .0001;
niter = scalar maximum number of iterations; default = 20;
options = matrix of n rows (same number of rows as pin) containing
column 1: desired fractional precision in parameter estimates.
Iterations are terminated if change in parameter vector (chg) on two
consecutive iterations is less than their corresponding elements
in options(:,1).  [ie. all(abs(chg*current parm est) < options(:,1))
on two consecutive iterations.], default = zeros().
column 2: maximum fractional step change in parameter vector.
Fractional change in elements of parameter vector is constrained to be
at most options(:,2) between sucessive iterations.
[ie. abs(chg(i))=abs(min([chg(i) options(i,2)*current param estimate])).],
default = Inf*ones().

OUTPUT VARIABLES
f = column vector of values computed: f = F(x,p).
p = column vector trial or final parameters. i.e, the solution.
kvg = scalar: = 1 if convergence, = 0 otherwise.
iter = scalar number of iterations used.
corp = correlation matrix for parameters.
covp = covariance matrix of the parameters.
covr = diag(covariance matrix of the residuals).
stdresid = standardized residuals.
Z = matrix that defines confidence region (see comments in the source).
r2 = coefficient of multiple determination.

All Zero guesses not acceptable