function dampedharmonicoscillator

tspan = [ 0, 10 ];
y0    = [ 0, 1 ]; %this is [x(0),x'(0)]
opts = odeset('Refine', 8);
[ tplot, yplot ] = ode45( @odefun, tspan, y0, opts );
figure(1), plot( tplot, yplot(:,1), '-k', ...
                 tplot, yplot(:,2), '--k' )
set(gca,'fontsize',16);
%phase portrait
figure(2), plot( yplot(:,1), yplot(:,2), '-k' )
set(gca,'fontsize',16);

function f = odefun(t,y)
omega = 2*pi;
zeta  = 0.1;  % parameters hard-coded
f = [ 0, 1; -omega^2, -2*zeta*omega ]*y;
end
end

