function ShowLSQNONLIN()
% Exponential fitting with LSQNONLIN

% Generate data (t(1),y(1)),...,(t(m),y(m)) where
% y(i) = a(1)*exp(-lambda(1)*t(i)) +...+a(p)*exp(-lambda(p)*t(i)) + noise

p = 2; % Number of exponentials
n = 2*p; % Number of Unknowns
a = [2;3]; % The coefficients
lambda = [0.5 ; 1.5]; % The time constants
m = 50; % The number of data points
tL = 0; % The lower bound of the sampling range
tU = 10; % The upper bound of the sampling range
t = linspace(tU,tL,m)'; % The sample times
tau = .1; % Noise factor
Noise = tau*randn(m,1); % Noise vector
y = exp(-t*lambda')*a + Noise; % The Data vector

% Generate the least square fit...
x0 = [1 1 1 2]'; % Starting vector
LB = [0 0 0 0]'; % Lower bound vector
UB = [inf inf inf inf ]'; % Upper bound vector
LB = []; UB = []; % This choice means no constraints imposed
options = optimset('Algorithm','levenberg-marquardt');
% type "optimset('lsqnonlin')" to see the full range of possibilities
% The call....
[xStar,resNorm,resVec,EXITFLAG,OUTPUT] = ...
lsqnonlin(@(x) MyR(x,t,y),x0,LB,UB,options);
% Extrtact the optimal a-values and lambda-values...
aStar = xStar(1:p);
lambdaStar = xStar(p+1:n);

% Display the data...
plot(t,y,'ro')
% and the fitting function...
N = 200;
tVals = linspace(0,tU,N)';
hold on
plot(tVals,exp(-tVals*lambdaStar')*aStar,'k','linewidth',2)
hold off
title(sprintf('a(1) = %5.3f lambda(1) = %5.3f a(2) = %5.3f lambda(2) = %5.3f',...
aStar(1),lambdaStar(1),aStar(2),lambdaStar(2)),'fontsize',14)
xlabel(sprintf('norm(R(xStar)/m = %5.3f MyR-evals = %3d',sqrt(resNorm)/m,OUTPUT.funcCount),'fontsize',14)
shg

function R = MyR(x,t,y)
% t and y are column m-vectors
% x is a column n-vector where n is even
% R is a column m-vector where
% R(i) = y(i) - (a(1)exp(-lambda(1)*t(i)) +...+a(p)exp(-lambda(p)*t(i)))
% with a = x(1:p) and lambda = x(p+1:n) where p = n/2.
% (The objective function if r(x) = .5*R'*R but lsqnonlin wants R)
n = length(x);
a = x(1:n/2);
lambda = x(n/2+1:n);
R = y - exp(-t*lambda')*a;