using System;
using static System.Math;
using static System.Console;
using System.Collections;
using System.Collections.Generic;

public partial class ode{

public static (vector,vector)
steprk23( Func<double,vector,vector> F, double x, vector y, double h)
{
	var k0 = F(x,y);
	var k1 = F(x+h/2  , y+(h/2  )*k0);
	var k2 = F(x+3*h/4, y+(3*h/4)*k1);
	var ka = (2*k0+3*k1+4*k2)/9;
	var kb = k1;
	var yh = y+ka*h;
	var er = (ka-kb)*h;
	return (yh,er);
}

public static (vector,vector) step2
(Func<double,vector,vector> F, double xo, vector yo, double xi, vector yi, double h)
{
	var ki=F(xi,yi);
	var c=(yo-yi+ki*(xi-xo))/(xi-xo).pow(2);
	var step=h/3;
	var t=xi+step;
	var yt=yi+ki*step+c*(step*step);
	var f2=F(t,yt);
	var d=(f2-ki-(2*(t-xi))*c)/(2*(t-xi)*(t-xo)+(t-xi).pow(2));
	var dy=d*(h*h*(xi+h-xo));
	var y=yi+ki*h+c*h*h+dy;
	return (y,dy/2);
}

public static vector driver
(
	Func<double,vector,vector> F, /* equation */
	double a, vector ya, /* initial condition: {a,y(a)} */
	double b, 
	double acc=1e-3,
	double eps=1e-3,
	double h=0.1, 
	List<double> xlist=null,
	List<vector> ylist=null,
	int limit=999
){
int nsteps=0;
if(xlist!=null) {xlist.Clear(); xlist.Add(a);}
if(ylist!=null) {ylist.Clear(); ylist.Add(ya);}
vector yh,dy,yo=ya.copy();
double xo=a;
do{
	if(a>=b) return ya;
	if(nsteps>limit) {
		Error.Write($"ode.driver: nsteps>{limit}\n");
		return ya;
		}
	if(a+h>b) h=b-a;
	if(nsteps==0) (yh,dy) = steprk23(F,a,ya,h);
	else (yh,dy) = step2(F,xo,yo,a,ya,h);
	double tol=(acc+yh.norm()*eps)*Sqrt(h/(b-a));
	double err=dy.norm();
	if(err<tol){
		nsteps++;
		xo=a; yo=ya.copy();
		a+=h; ya=yh;
		if(xlist!=null) xlist.Add(a);
		if(ylist!=null) ylist.Add(ya);
		}
	else { Error.WriteLine($"driver: bad step at {a}"); }
	if(err>0) h = h*Pow(tol/err,0.25)*0.95;
	else h=2*h;
	}while(true);
}// driver
}// class