tspan=[0,200];       % time span
x0 = [2,2];          % initial position and velocity
[t,x] = ode45(@harmonicoscillator,tspan, x0); 
figure(1)
hold off 
plot(x(:,1),x(:,2)); % phase portrait
xlabel('position','FontSize',10)
ylabel('velocity','FontSize',10)
title('Forced and dumped harmonic oscillator','FontSize',10);
figure(2) 
hold off  
comet(x)            % trajectory
xlabel('time','FontSize',10)
ylabel('position','FontSize',10)
title('Forced and dumped harmonic oscillator','FontSize',10); 
set(figure(3),'Position', [500 100 400 600]);
clf;
subplot(2,1,1)
plot(x(:,1),x(:,2)); % phase portrait
xlabel('position','FontSize',10)
ylabel('velocity','FontSize',10)
subplot(2,1,2)
comet(x)            % trajectory
xlabel('time','FontSize',10)
ylabel('position','FontSize',10)
