%% Model Population Data
% Use population data at 1910, 1950, and 1990 to model population with function P(t)=a_1+a_2e^{a_3t}.
% In part a newtonsys does not converge. The numbers are too large. So, one rescales the data. Use t=(year-1900)/100.
%                     \matrix{\mbox{Years} & .1 & .3 & .5 & .7 & .9 \cr \mbox{Pop.} & .92 & 1.228 & 1.507 & 2.050 & 2.487}

x1 = [1;1;1];   % column vector!
x = newtonsys(@nlsystem,x1);

%
% Display converged solution, plot model function with data, and find P(2000).

a=x(:,end)
P=@(t) 100*(a(1)+a(2)*exp(a(3)*(t-1900)/100));
P(1950)

x2=[1910, 1930, 1950, 1970, 1990];
y=[92, 122.8, 150.7, 205, 248.7];
plot(x2,y,'o')
xlabel('year')
ylabel('Population')
hold on
fplot(P,[1910,1990])
hold off
P(2000)

%
%

function [f,J] = nlsystem(a)
    f = zeros(3,1);   % ensure a column vector output
    f(1) = a(1)+a(2)*exp(a(3)*.10)-.92;
    f(2) = a(1)+a(2)*exp(a(3)*.50)-1.507;
    f(3) = a(1)+a(2)*exp(a(3)*.90)-2.487;
    J(1,:) = [1, exp(a(3)*.10), .10*a(2)*exp(a(3)*.10)];
    J(2,:) = [1, exp(a(3)*.50), .50*a(2)*exp(a(3)*.50)];
    J(3,:) = [1, exp(a(3)*.90), .90*a(2)*exp(a(3)*.90)];
end

function x = newtonsys(f,x1)
% NEWTONSYS   Newton's method for a system of equations.
% Input:
%   f        function that computes residual and Jacobian matrix
%   x1       initial root approximation (n-vector)
% Output       
%   x        array of approximations (one per column, last is best)

% Operating parameters.
funtol = 1000*eps;  xtol = 1000*eps;  maxiter = 40;

x = x1(:);  
[y,J] = f(x1);
dx = Inf;
k = 1;

while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
  dx = -(J\y);   % Newton step
  x(:,k+1) = x(:,k) + dx;

  k = k+1;
  [y,J] = f(x(:,k));
end

if k==maxiter, warning('Maximum number of iterations reached.')
end
end