clear all;
global m;
global b;
global k;
m=1;
b=0.5;
k=1;
 
 
function y=ujfun(x)
global m;
global b;
global k;
   y(1)= x(2);
   y(2)=-k/m*x(1)-b/m*x(2);
endfunction;
 
x0=[1;0];
t=[0:0.01:20];
Z=lsode("ujfun", x0, t);
 
w=sqrt(k/m-(b/(2*m))^2)
subplot(1,2,1)
plot(t,Z(:,1),'r*',t,exp(-b/(2*m)*t).*(cos(w*t)+1/2*b/(m*w)*sin(w*t)),'b')
subplot(1,2,2)
 
% portret fazowy
[X,Y]=meshgrid (-1:0.1:1, -1:0.1:1);
F=Y;
G=-k/m*X-b/m*Y;
quiver(X,Y,F,G)
hold on 
plot(Z(:,1),Z(:,2),'r')