% Integrate a spring-mass system

k=1; beta=.01;

X0 = [1;3;k;beta];
T = linspace(0,1000,100000);
sol = lsode("forced_pendulum", X0, T);

% Time series plot
%plot(sol(:,1:2));

% Phase diagram - velocity as function of position
title "The phase space diagram of a damped pendulum"
xlabel "Position"
ylabel "Velocity"
% Skip the initial segment of the trajectory
offset=70000;
len = rows(sol);
plot(sol(offset:len,1),sol(offset:len,2));
