function ShowMarkov()
% Looks at dy/dt = Q*y(t) where Q is an ``intensity matrix''.
% This is a matrix whose column sums are zero and whose off-diagonal
% entries are positive.

% Generate a random intensity matrix...
n = 8;
Q = rand(n,n);
for k=1:n
Q(k,k) = 0;
Q(k,k) = -sum(Q(:,k));
end

clc
disp('Look at exp(Q/3)...')
disp(' ' )
disp(expm(Q/3))
% Generate a random initial condition for dy/dt = Q*y(t)...
y0 = rand(n,1); y0 = y0/sum(y0);
% The solution is given by y(t) = exp(Qt)*y0, The matrix
%
% exp(Qt) = I + Qt + (Qt)^2/2! + (Qt)^3/3! +...
%
% is a probability matrix. (Nonnegative entries and column sums equal to
% one.)

% Display the solutions over the interval [0,1]
m = 401;
tVals = linspace(0,4,m);
[t,y] = ode45(@(t,y) MyF(t,y,Q),tVals,y0);
plot(t,y)
shg

% Reverse integrate for various t-values to see if we can recover the
% initial condition. Notice the way the interval of integration is
% specified.

disp(' ')
disp('Starting time Error in recovering y0')
for k=[401 301 201 101 51]
[t,yRev] = ode45(@(t,y) MyF(t,y,Q),[tVals(k) 0],y(k,:)');
err = norm(y0-yRev(length(t),:)');
fprintf('%6.3f %10.3e\n',tVals(k),err)
end

function yp = MyF(t,y,Q)
yp = Q*y;