#include <kv/allsol.hpp>
#include <kv/constants.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct FiveSol {
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
		int i, j;
		ub::vector<T> y(4);
		ub::vector<T> f(4);
		ub::matrix<T> mT(4,4);
		ub::matrix<T> G(4,4);
		ub::vector<T> J(4);

		const static T rb = kv::constants<T>::str("10000");
		const static T rc = kv::constants<T>::str("5000");
		const static T vcc = kv::constants<T>::str("-5");
		const static T af = kv::constants<T>::str("0.99");
		const static T ar = kv::constants<T>::str("0.5");
		const static T is = kv::constants<T>::str("1e-9");
		const static T vt = kv::constants<T>::str("0.053");
		const static T vs = kv::constants<T>::str("-0.64");

		T gb, gc;

		gb = 1. / rb;
		gc = 1. / rc;

		for (i=0; i<4; i++) {
			for (j=0; j<4; j++) {
				mT(i,j) = 0.;
			}
		}
		mT(0,0) = 1.; mT(0,1) = -ar;
		mT(1,0) = -af; mT(1,1) = 1.;
		mT(2,2) = 1.; mT(2,3) = -ar;
		mT(3,2) = -af, mT(3,3) = 1.;

		G(0,0) = 2. * gb + gc;
		G(0,1) = -(gb + gc);
		G(0,2) = -2. * gb;
		G(0,3) = gb;
		G(1,0) = -(gb + gc);
		G(1,1) = gb + gc;
		G(1,2) = gb;
		G(1,3) = 0.;
		G(2,0) = -2. * gb;
		G(2,1) = gb;
		G(2,2) = 2. * gb + gc;
		G(2,3) = -(gb + gc);
		G(3,0) = gb;
		G(3,1) = 0.;
		G(3,2) = -(gb + gc);
		G(3,3) = gb + gc;

		J(0) = gc * vcc;
		J(1) = gb * vs - gc * vcc;
		J(2) = gc * vcc;
		J(3) = gb * vs - gc * vcc;

		f(0) = is / af * (exp(x(0) / vt) - 1.);
		f(1) = is / ar * (exp(x(1) / vt) - 1.);
		f(2) = is / af * (exp(x(2) / vt) - 1.);
		f(3) = is / ar * (exp(x(3) / vt) - 1.);

		y = prod(mT, f) + prod(G, x) + J;

		return y;
	}
};


int main()
{
	int i;
	ub::vector<itv> I;

	std::cout.precision(17);

	I.resize(4);
	for (i=0; i<I.size(); i++) I(i) = itv(-10., 10.);

	kv::allsol(FiveSol(), I);
}