function rk
tspan = [0;50];
y0 = [6;0];

options = odeset('RelTol',1e-5,'AbsTol',1e-4);
[t,y] = ode45(@fu1,tspan,y0,options);

plot(t,y(:,1),'.-');
figure
plot(y(:,1),y(:,2),'.-');

% --------------------------------------------------------------------------
function ydot = fu1(t,y)

w = 3;
b=.5;
ydot = [y(2); -b*y(2)-w^2*y(1)];