using System;
using static System.Console;
using static System.Math;
using System.Collections.Generic;
public static partial class ode{

public static vector driver(// ode driver
	Func<double,vector,vector> F, /* equation */
	double a, vector ya, double b, 
	double acc=0.01, double eps=0.01, double h=0, 
	List<double> xlist=null, List<vector> ylist=null,
	int limit=999
	){
int nsteps=0;
if(h==0)h=(b-a)/64;
if(xlist!=null) {xlist.Clear(); xlist.Add(a);}
if(ylist!=null) {ylist.Clear(); ylist.Add(ya);}
double x=a;
var y=ya.copy();
int n=y.size;
var tol=new vector(n);
do{
	if(x>=b) return y;
	if(nsteps>limit) throw new Exception($"driver: nsteps>{limit}\n");
	if(x+h>b) h=b-x;
	var (yh,erv) = rkstep23(F,x,y,h);

	for(int i=0;i<n;i++) tol[i]=Max(acc,Abs(yh[i])*eps)*Sqrt(h/(b-a));
	bool ok=true; for(int i=0;i<n;i++) if(Abs(erv[i])>tol[i])ok=false;

	if(ok){
		x+=h; y=yh; nsteps++;
		if(xlist!=null) xlist.Add(x);
		if(ylist!=null) ylist.Add(y);
		}

	double factor = tol[0]/Abs(erv[0]);
	for(int i=1;i<n;i++) factor=Min(factor,tol[i]/Abs(erv[i]));
	h *= Min( Pow(factor,0.25)*0.95 ,2);

	}while(true);

}//driver

}//class