mode(-1);

getf('FixedRK.sci');
getf('RKstep.sci');
getf('funcs.sci');
%ODEOPTIONS=[2,0,0,%inf,0,2,500,12,5,0,-1,-1];

ok=1;
while ok do
  str=x_mdialog(['Résolution d''un problème à trois corps'; ...
	  'Comparaison de méthodes à pas fixes et d''une méthode avec contrôle du pas'], ...
      ['Ordre (0 pour variable)'; 'Nombre de pas de temps'; 'Temps final'], ...
      ['4'; '5000'; '17.1']);

  ok = (str~=[]);
  if ok then
    ordre=evstr(str(1)); 
    n=evstr(str(2));
    tf=evstr(str(3));
    t0=0;
    y0=[0.994; 0; 0; -2.0015851063791];
    h=tf/n;
  
    if ordre ~=0 then
      [tvals, yvals] = FixedRK('aren', t0, y0, h, ordre, n);
      yvref= ode(y0, t0, tf, aren);
      xbasc();
      plot2d(yvals(:,1), yvals(:,2), style=1);
      plot2d(yvref(2,:), yvref(3,:), style=2, leg='Pas variable');
      xtitle('Solutions du problème à trois corps');
    else
      yvref= ode(y0, t0, tf, aren);      
      nts= size(yvref,2); ts=linspace(t0, tf, nts)';
      x_message(['évolution du pas de temps et de la solution'; string(nts)+' pas de temps']);    
      xbasc();
      dts=yvref(1,2:nts)-yvref(1,1:nts-1);
      subplot(311); plot2d(ts(2:$), dts);
      subplot(312); plot2d(ts, yvref(2,:));
      subplot(313); plot2d(ts, yvref(3,:));
    end
  end
end