#include <iostream>
#include <kv/ode.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
	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) = -x(0); 
		return y;
	}
};

int main()
{
	ub::vector<itv> x;
	itv end;
	int r, i;
	ub::vector< kv::psa<itv> > result_psa;

	std::cout.precision(17);

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

	end = 10.;
	r = kv::ode(Func(), x, itv(0.), end, kv::ode_param<double>().set_order(5), &result_psa);

	if (r == 0) {
		std::cout << "Cannot calculate solution.\n";
	} else if (r == 1) {
		std::cout << "Solution calculated until t = " << end << ".\n";
		std::cout << x << "\n";
	} else {
		std::cout << "Solution calculated.\n";
		std::cout << x << "\n";
	}

	for (i=0; i<result_psa.size(); i++) {
		std::cout << result_psa(i) << "\n";
	}
}