% 4-15-2008
% Jacobi redux

clear all;
close all;

% Parameters for the iteration
max_iterations = 1000;
TOLERANCE = 10^(-6);

% Initial conditions, note that here we are telling matlab that we want a
% matrix rather than a vector
x(1,1:41) = 1;
y(1,1:41) = 2;
z(1,1:41) = 2;

% Vector of parameter values we would like to use
c = 6:0.1:10;

% This will keep track of what the min/max number of iterations that our
% solver has used over all of the parameter values c
i_min = 100;
i_max = 0;

% Parameter loop
for j=1:length(c)
    % Linear system solver loop
    for i=2:max_iterations
        % Jacobi update
        x(i,j) = (c(j) + y(i-1,j) - z(i-1,j)) / 4;
        y(i,j) = (21 + 4*x(i-1,j) + z(i-1,j)) / 8;
        z(i,j) = (15 + 2*x(i-1,j) - y(i-1,j)) / 5;

        % Check to see if we have reached convergence
        if sqrt( (x(i,j) - x(i-1,j))^2 + (y(i,j) - y(i-1,j))^2 ...
                + (z(i,j) - z(i-1,j))^2 ) < TOLERANCE
            break;
        % Now check to see if our solver may have failed, i.e. one of our
        % values has blown up.  Note the use of the logical 'or' operator
        elseif abs(x(i,j)) > 10^(12) || abs(y(i,j)) > 10^(12) || abs(z(i,j)) > 10^(12)
            break;
        end
    end
    
    % Check to see if we used less iterations to solve the loop this time
    % than all the past times
    i_min = min(i_min,i);
    % Same as above but for the maximum
    i_max = max(i_max,i);
end