function [t,y] = linsys(A,x0,T,siz)
% Input: Matriisi A, alkuarvovektori x0, loppuaika T (aikaväli: 0 .. T)
% Valinnainen siz: kuinka moneen osaan aikaväli jaetaan, oletus 100.
%
% Output: t - vektori:  diskretoitu aika-akseli (oletus 100-pituinen)
%         y - matriisi: 100 x n, kukin sarake edustaa ratkaisufunktion
%                       arvoja t-aikapisteissä.
% 
% Esim: A=[1 0 0;1 3 0;0 1 1];x0=[1;2;3];
%       [t,y]=linsys(A,x0,2);
%       plot(t,y)                       % aikariippuvuusparvi (mieliv. n)
%       plot(y(:,1),y(:,2))             % faasitaso (1,2)  
%       plot3(y(:,1),y(:,2),y(:,3))     % faasiavaruus (3-d)
%       plot(y(:,i),y(:,j))             % projektio (i,j)-faasitasossa
%                                       % (jos n > 2)
if nargin < 4, siz = 100; end; 
t = linspace(0,T,siz);
m=size(A);m=m(1);
y = zeros(m,siz);
for i=1:siz,
y(:,i) = expm(A * t(i)) * x0;
end;
y=y';      % Transponoidaan y-matriisi,jotta yhdenmukainen ode-funkt. kanssa
