#include <cmath>
#include <iostream>
#include <tuple>
#include <vector>

#include "euler_heun.h"
#include "maple.h"

/********************************
 * MUST BE COMPILED WITH C++ 11 *
 ********************************/

// Funciton declarations
int main();
double f0( double t, double y );
double f1( double t, double y );

// Funciton definitions
double f0( double t, double y ) {
    return (t - 1)*y + 0.5;
}

double f1( double t, double y ) {
    double const PI{ std::acos(-1.0) };
    return -y + 10.0/std::sqrt( PI )*std::exp( -100.0*(t - 2.0)*(t - 2.0) );
}

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

    std::cout << "# Adaptive Euler-Heun" << std::endl;
    
    std::cout << std::endl << "# y'(t) = (t - 1)y(t) + 0.5" << std::endl;
    std::cout << "#  y(0) = 1.2" << std::endl;
    std::cout << "#  y(2) = ?" << std::endl << std::endl;
    std::cout << "# Maximum error of 0.5:" << std::endl;

    auto result = adaptive_euler_heun( f0, 0, 1.2, 2, 0.01, std::make_pair( 1e-10, 0.5 ), 0.5/2/2, 10000 );

    vector_size_maple( std::get<0>( result ) );
    vector_print_maple( " t0", std::get<0>( result ) );
    vector_print_maple( " y0", std::get<1>( result ) );
    vector_print_maple( "dy0", std::get<2>( result ) );

    std::cout << std::endl << "# Maximum error of 0.1:" << std::endl;

    result = adaptive_euler_heun( f0, 0, 1.2, 2, 0.01, std::make_pair( 1e-10, 0.5 ), 0.25, 10000 );

    vector_size_maple( std::get<0>( result ) );
    vector_print_maple( " t1", std::get<0>( result ) );
    vector_print_maple( " y1", std::get<1>( result ) );
    vector_print_maple( "dy1", std::get<2>( result ) );

    std::cout << std::endl << "#                                  2" << std::endl;
    std::cout << "#                 10   -100 (t - 2)" << std::endl;
    std::cout << "# y'(t) + y(t) = ---- e" << std::endl;
    std::cout << "#                \\/pi" << std::endl;
    std::cout << "#  y(0) = 1.2" << std::endl;
    std::cout << "#  y(2) = ?" << std::endl << std::endl;
    std::cout << "# Maximum error of 0.5:" << std::endl;

    result = adaptive_euler_heun( f1, 0, 1.2, 4, 0.01, std::make_pair( 1e-10, 0.25 ), 1.0, 10000 );

    vector_size_maple( std::get<0>( result ) );
    vector_print_maple( " t2", std::get<0>( result ) );
    vector_print_maple( " y2", std::get<1>( result ) );
    vector_print_maple( "dy2", std::get<2>( result ) );

    std::cout << std::endl << "# Maximum error of 0.1:" << std::endl;

    result = adaptive_euler_heun( f1, 0, 1.2, 2, 0.01, std::make_pair( 1e-10, 0.5 ), 0.1, 10000 );

    vector_size_maple( std::get<0>( result ) );
    vector_print_maple( " t2", std::get<0>( result ) );
    vector_print_maple( " y2", std::get<1>( result ) );
    vector_print_maple( "dy2", std::get<2>( result ) );

    return 0;
}