import org.apache.commons.math3.ode.FirstOrderDifferentialEquations
import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator
import org.apache.commons.math3.ode.ContinuousOutputModel

class vanderpol implements FirstOrderDifferentialEquations{

	Closure f;
	double mu,om;

	vanderpol(double mu, double om, Closure f) {
		this.f = f;
		this.om = om;
		this.mu = mu;
 	}

	int getDimension() {
		return 2;
	}

	void computeDerivatives(double t, double[] y, double[] yDot) {
		yDot[0] = y[1]
		yDot[1] = mu*(1-y[0]*y[0])*y[1]-om*om*y[0]+f(t)
	}

}

double dt=0.01
CRK = new ClassicalRungeKuttaIntegrator(dt)

f={0.0}

ode = new vanderpol(1.0,1.0,f);
double[] y = [0.0, -0.01]; // initial state

tracker = new ContinuousOutputModel()

CRK.addStepHandler(tracker);

double tot=150.0
CRK.integrate(ode, 0.0, y, tot, y);
u=[]
v=[]
(0..tot/dt).each{
	tracker.setInterpolatedTime(it*dt)
	res=tracker.getInterpolatedState()
	u.add(res[0])
	v.add(res[1])
	//println it*dt+" "+res[0]+" "+res[1]
}

thePlot.clear()
//thePlot.addFunction(new plotfunction(dt,u as double[]))
thePlot.addFunction(new plotfunction(u as double[],v as double[]))
thePlot.setAutoColor(true)
thePlot.xLabel("u")
thePlot.yLabel("v")
thePlot.Title("van der Pol oscillator")
thePlot.setMarker(false)
thePlot.show()