#include <kv/poincaremap.hpp>
#include <kv/kraw-approx.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct VDP {
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
		ub::vector<T> y(2);

		y(0) = x(1);
		y(1) = 1. * (1. - x(0)*x(0)) * x(1) - x(0);

		return y;
	}
};

struct VDPPoincareSection {
	template <class T> T operator() (const ub::vector<T>& x){
		T y;

		y = x(0) - 0.;

		return y;
	}
};


int main()
{
	ub::vector<double> x;
	ub::vector<itv> ix;
	bool r;

	std::cout.precision(17);

	VDP f;
	VDPPoincareSection ps;
	kv::PoincareMap<VDP,VDPPoincareSection,double> po(f, ps, (itv)0.);

	x.resize(3);
	x(0) = 0.; 
	x(1) = 2.;
	x(2) = 6.;

	r = kv::krawczyk_approx(po, x, ix, 5);
	if (r) {
		std::cout << "solution found.\n";
		std::cout << ix << "\n";
	}
}