function F = MyExpM(T)
% T is a normalized upper quasi-triangular matrix with distinct eigenvalues.
% F = exp(T)

[n,n] = size(T);
if n==1
% 1-by-1 case...
F = exp(T);
elseif n==2
% 2-by-2 case
if T(2,1)~=0
% Complex eigenvalues...
F = exp2(T);
else
% Real eigenvalues...
F11 = exp(T(1,1));
F22 = exp(T(2,2));
F12 = T(1,2)*(F22 - F11)/(T(2,2) - T(1,1));
F = [F11 F12; 0 F22];
end
else
% Write
% T11 T12 F11 F12
% T = F =
% 0 T22 0 F22
%
% where T11 is either 1-by-1 or 2-by-2. Regardless, we compute
% F11 = exp(T11), F22 = exp(T22) (recursively), and then solve
% T11*F12 - F12*T22 = F11*T12 - T12*F22
% for F12...
if T(2,1)==0
% T11 is 1-by-1...
T11 = T(1,1); T12 = T(1,2:n); T22 = T(2:n,2:n);
F11 = exp(T11);
F22 = MyExpM(T22);
F12 = TinySylvesterSolve(T11,T22,F11*T12-T12*F22);
F = [F11 F12; zeros(n-1,1) F22];
else
% T11 is 2-by-2 ...
T11 = T(1:2,1:2); T12 = T(1:2,3:n); T22 = T(3:n,3:n);
F11 = exp2(T11);
F22 = MyExpM(T22);
F12 = TinySylvesterSolve(T11,T22,F11*T12-T12*F22);
F = [F11 F12; zeros(n-2,2) F22];
end
end

function F = exp2(T)
% F = exp(T) where T = [lambda alpha; beta lambda]
% and alpha*beta<0
%

lambda = T(1,1);
alfa = T(1,2);
beta = T(2,1);
% Since F = exp(lambda)*exp([0 alpha ; beta 0]) we need an expression
% for exp([0 alpha ; beta 0]), i.e., the solution to

% Xdot = [0 alfa; beta 0] X X(0) = I
% Comparing entries we see
% X11_dot = alfa*X21
% X21_dot = beta*X11
% from which we conclude that
% X11_dotdot = -mu*X11 X11(0)=1, X11_dot(0)=0
% where
% mu = sqrt(-alfa*beta).
% This has solution
% X11(t) = cos(mu*t)
% Likewise
% X21_dotdot = -tau*X21 X21(0) = X21_dot(0)=beta
% has solution
% X21(t) = sign(beta)sqrt(-beta/alfa)*sin(mu*t)
% Similar reasoning reveals
% X22(t) = cos(mu*t)
% and
% X12(t) = sign(alfa)sqrt(-alfa/beta)*sin(mu*t)

mu = sqrt(-alfa*beta); c = cos(mu); s = sin(mu);
F = exp(lambda)*[c sign(alfa)*sqrt(-alfa/beta)*s;...
sign(beta)*sqrt(-beta/alfa)*s c];