function plotsolution(A,x0,tmax,dt,range)
% plots a solution for the 2-D linear system defined by A and the initial
% condition *column* vector x0. (it has to be a column vector
% because of the use of 'inv' below.) 
% tmax is maximum time and dt is time step used for plotting.
% range is the maximum x and y range for plotting.

% test for size
[m,n] = size(A);
if m ~= 2 | n ~= 2
   'Matrix is not 2-D - quitting.'
   return
end

% eigenanalysis and test for real result
[c,l] = eig(A);
Ca = c(:,1);
Cb = c(:,2);
La = l(1,1);
Lb = l(2,2);
if ~isreal(c)
   'Eigenvectors/values are imaginary - quitting.'
   return
end

% find the values of alpha corresponding to initial condition x0
alpha = inv(c)*x0;

% plot solution for time range 0:tmax in steps of dt
t = 0:dt:tmax;
x = alpha(1)*Ca*exp(La*t) + alpha(2)*Cb*exp(Lb*t);

% truncate if needed
imax = length(t);
for i = 1:length(t),
   if x(1,i) > range | x(1,i) < -range
      imax = i;
      break;
   elseif x(2,i) > range | x(2,i) < -range
      imax = i;
      break;
   end
end
      
hold on;
plot(x(1,1:imax),x(2,1:imax),'m-');