#include <iostream>
#include "shooting.h"
#include "vec.h"

// Function declarations
vec<2> f( double x, vec<2> u );
int main();

// Function definitions

// A linear ODE:
//     u"(x) - u'(x) - 3.2u(x) = -1
vec<2> f1( double x, vec<2> u ) {
    return vec<2>{u[1], 3.2*u[0] + u[1] - 1.0};
}

// A non-linear ODE:
//     u"(x) - u'(x)(u(x) - 3) = 0
vec<2> f2( double x, vec<2> u ) {
    return vec<2>{u[1], u[1]*(u[0] - 3.0)};
}

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

    unsigned int const N{ 16 };
      ///////////////
     // Example 1 //
    ///////////////

    double   a{0.0};
    double   b{1.0};
    double u_a{4.7};
    double u_b{3.8};

    ivp<vec<2>> soln1{ bvp( f1,
                            std::make_pair(   a,   b ),
                            std::make_pair( u_a, u_b ),
                            1e-10, 10, 1e-2, std::make_pair( 1e-5, 1e-1 ) ) };

    double h{(1.0 - 0.0)/N};
    std::cout << "xs = [" << a;

    for ( unsigned int k{1}; k < N; ++k ) {
        std::cout << ", " << (a + k*h);
    }

    std::cout << ", " << b << "];" << std::endl;

    std::cout << "us = [" << soln1( a )[0];

    for ( unsigned int k{1}; k <= N; ++k ) {
        std::cout << ", " << soln1( a + k*h )[0];
    }

    std::cout << "];" << std::endl;
 
    std::cout << "plot( xs, us );" << std::endl;

      ///////////////
     // Example 2 //
    ///////////////

    // The exact solution is
    //  u(x) = 3-1.543404638418209*tanh(0.7717023192091042*x-1.543404638418208)

    a = 1.0;
    b = 2.0;
    u_a = 4.0;
    u_b = 3.0;

    ivp<vec<2>> soln2{ bvp( f2,
                            std::make_pair(   a,   b ),
                            std::make_pair( u_a, u_b ),
                            1e-10, 20, 1e-2, std::make_pair( 1e-5, 1e-1 ) ) };

    h = (b - a)/N;
    std::cout << "xs = [" << a;

    for ( unsigned int k{1}; k < N; ++k ) {
        std::cout << ", " << (a + k*h);
    }

    std::cout << ", " << b << "];" << std::endl;

    std::cout << "us = [" << soln2( a )[0];

    for ( unsigned int k{1}; k <= N; ++k ) {
        std::cout << ", " << soln2( a + k*h )[0];
    }

    std::cout << "];" << std::endl;
 
    std::cout << "plot( xs, us );" << std::endl;

    return 0;
}