% MalateEqs Script solves the nonlinear equations in the malate problem.
% Additional M-files required:
%   X0.m, eqlbrm.m, eqljac.m
% 
format short e;
format compact;
% 
%	Next, constants and data
% 
epsk=1.0/sqrt(10^(4.05)); 
c = epsk*10^(2.24);
% epsk = 0.00944; c = 1.64;
% 
% Always a(2) = 156.
% 
a = [0;156.0;0];
% 
%	Values of a(3) for 3 experiments
% 
a3=[100 200 400];
% 
%  Values of a(1) (set), 3 experiments
% 
a1 =     [31.25   62.5   125.0    187.5    250.0    500.0    ]';
a1 = [a1 [15.0    31.25   62.5    125.0    250.0    500.0    ]'];
a1 = [a1 [15.0    31.25   62.5    125.0    250.0    500.0    ]'];
% 
% LOOP: for each experiment in three series of six,
% solve the equilibrium equations for {x,y,z}.
% 
xtol = 1e-6; rtol = 1e-6; imax = 10;
xsol = zeros(3,6);
for m=1:3
  a(3)=a3(m);
  disp(sprintf('\n Run %d of 6 experiments',m))  % Table Header
  disp('|    x     |    y     |    z     |')
  disp('----------------------------------')
  for k=1:6
    a(1)=a1(k,m);
%   Newton's method: set initial guess, loop control.
    x = X0(a,epsk,c);
    snorm = 1 + xtol; fnorm = 1 + rtol; i = 0;
    while (snorm > xtol & fnorm > rtol & i <= imax)
      r = eqlbrm(x,a,epsk,c);   % Residual
      J = eqljac(x,a,epsk,c);   % Jacobian
      s = - (J \ r);            % Newton correction
      x = x + s;                % Update root approx, loop controls
      snorm = norm(s);
      rnorm = norm(r);
      i = i + 1;
    end
% 
% Display the results.
% 
    disp(sprintf(' %10.4e %10.4e %10.4e',x))
    xsol(:,k) = x;
  end
end