function  DglDemo(A,np,nt,r,varargin)
% DglDemo(A,np,nt,r)  shows solution curves of linear ODEs
% size(A) = [2,2];  np = number of points;  nt = time steps
% TRY:  DglDemo;  DglDemo([0 -1; 1 0]);  DglDemo([-0.1 2; 2 -0.4]);

switch nargin     %.. Defaults for varying number of inputs
 case 0
np=60; nt=30; A=[ -0.9  -1; 1 -0.8]; r=1;
 case 1
 np=60; nt=30; r=1;
 case 2
 nt=30; r=1;
 otherwise
end
% ..........Preparation of figure ................
figure, axis off, axis([-1  1  -1   1 ] ), axis equal,hold on
set(gcf,'Color','black');

% ..........Initial points and their iteration .......
rip = r*(rand(2,np) - 0.5 );       % np random initial points
pkt=rip; qkt=rip;
    plot(pkt(1,:),pkt(2,:),'o','erase','none')       
    drawnow
B=expm(A/nt);   C=expm(-A/nt); % Matrix-Exponential and Inverse
 for k=1:nt                    % nt time steps in this loop
    pkt=B*pkt; qkt=C*qkt; % Execute one time-step forward& backward
	plot(pkt(1,:),pkt(2,:),'.y','erase','none')   % forward
	plot(qkt(1,:),qkt(2,:),'.m','erase','none')   % backward 
    drawnow
 end
% =================== alternatives ==========================

% for k=1:nt
%   point(k)=line('LineStyle', 'none', 'marker' ,'.', ...
%       'color','y', 'EraseMode','none','xdata',[],'ydata',[]);
% end
%......%  set(point(k),  'xdata',pkt(1,:),  'ydata',pkt(2,:)),
