function BonusOctave

%size of your phase space
n=30;
%initial time
t_ini=0
%final time
t_fin=60
%frames to generate
frames=100
%some settings
set(0, 'defaultfigurevisible', 'off');

%your trajectories
x=zeros(frames,n*n*2);

%calculating trajectories for each entry in your phase space
it=1;
for i=0:1:n-1
  for j=0:1:n-1
    %calculate the trajectory of the initial condition (q_i,p_j)
    X=lsode("f",[-.1+.2/(n-1)*i;-.1+.2/(n-1)*j],(t=linspace(t_ini,t_fin,frames)'));
    x(:,it)=X(:,1);x(:,it+1)=X(:,2);
    it+=2;
  end
end

%Preparing for visualisation
y=zeros(frames,2);
graphics_toolkit("qt"); %a couple of years ago it also worked with "gnuplot", try it if "qt" does not work

%making the frames
for T=1:1:frames;
  for i =1:2:n*n*2
    y(i,1)=x(T,i);y(i,2)=x(T,i+1);
  end
  plot(y(:,1),y(:,2),".");
  axis([-3,3,-3,3]);
  filename=sprintf('output/%05d.png',T);
  print(filename, '-dpng');
end

end

%This is specified by your model
function xdot=f(x,t)
xdot(1)=x(2);
xdot(2)=-x(1)^3+x(1);
end


