function J=cstrJ(x)
% Reactor in series model, the Jacobian.
%
% x = [C(1),C(2),... C(N-1),beta]
% f(i) = beta C(i)^2 + C(i) - C(i-1)

global N CN

% Initialize matrix.
J =  zeros(N);

% First row.
J(1,1) = 2.0*x(N)*x(1)+1.0;
J(1,N) = x(1)^2.0;

% Middle rows.
for i=2:N-1
 J(i,i) = 2.0*x(N)*x(i)+1.0;
 J(i,i-1) = -1.0;
 J(i,N) = x(i)^2.0; 
end

% Last row.
J(N,N-1) = -1.0;
J(N,N) = CN^2.0;