#include <iostream>
#include <cmath>
#include "ivp.h"
#include "vec.h"

// Function declarations
int main();
vec<4> f( double x, vec<4> w );


// Function definitions
vec<4> f( double x, vec<4> w ) {
    return vec<4>{
        w[1],
        x + x*w[0],
        w[3],
        0 + x*w[2]
    };
}

int main() {
    std::cout.precision( 16 );
    std::clog.precision( 16 );

    ivp<vec<4>> u{ f, 0.0, vec<4>{ 1.2, 0.0, 0.0, 1.0 }, 0.1, 1e-10, 1e-5, vec<4>::norm };

    for ( double x{0.0}; x <= 1.0; x += 0.1 ) {
       std::cout << "u(" << x << ") = " << u(x) << std::endl;
    }

    return 0;
}