  function TestMyExp()
% Explores what m = 2^k should be so that R_{m}(A) is good enough for
% a given tol. 

% Key findings as a result of playing with this script:

% 1. n doesn't matter.
% 2. It makes sense to choose m = 2^k since high powers can
%    be reached cheaply.
% 3. The smallest tol that can be handled is 10^-10
clc
n = 100;
fprintf('What must k be so that R_{m}(A) is good enough where m = 2^k\n')
fprintf('\nn = %1d\n',n)
fprintf('\n  tol        k            m            Actual Error\n')
fprintf('------------------------------------------------------')
for tol = logspace(-1,-12,12)
   fprintf('\n%7.0e',tol)
   A = randIntensity(n);
   F0 = expm(A);
   k = 1;
   F = MyExpPrototype(A,1);
   itMax = 20;
   while max(max(abs(F - F0)))>tol && k<=itMax
      k=k+1;
      F = MyExpPrototype(A,k);
   end
   fprintf('  %6d  %12d %20.3e',k,2^k,max(max(abs(F - F0))))
   
end
fprintf('\n')

  function A = randIntensity(n)
% A is a random nxn normalized intensity matrix.
A = rand(n,n)/n;
for j=1:n
    A(j,j) = 0;
    A(j,j) = -sum(A(:,j));
end

  function F = MyExpPrototype(A,k)
% A is an nxn normalized intensity matrix with n<=1000
% F = R_{m}(A) where m = 2^k
[n,n] = size(A);
I = eye(n,n);
A0 = A/2^k;
F = (I -A0/2)\(I + A0/2);
% Repeated squaring...
for j=1:k
    F = F*F;
end

% 1 pt for evidence that your choice of m does not depend on n
% 1 pt for evidence that your choice of m works almost all the time
% Up to -1 for not testing the full range of tols and ns
% -1 for an error that makes results useless/incorrect
% Up to -1 for evidence not presented neatly as output (if code is understandable enough, no points lost)

