#include <math.h>

#include <iostream>
using namespace std; 
//
//#####################################################################
//
//                        EulerSample.cpp 
//
//#####################################################################
//
// This program solves the ordinary differential equation
// 
// dy/dt = f(y,t) y(t_0) = y_0
// 
// using Euler's method. Euler's method is implemented as a separate 
// routine where the function defining the ODE is passed in as a 
// function pointer. 
//
//
// Version 01/04/09
//#####################################################################
//
//
// ODE Specification
//
double f(double y, double t)
{
    return t*cos(y);
}

//
// Euler's method implementation declaration (prototype)
//
double EulerStep(double yn, double tn, double dt,
double(*F)(double,double));

int main()
{

    long    i;

    long    N;
    double dt;

    double yn; double tn;

    double tInitial =  0.0;
    double tFinal   =  2.0;
    double yInitial =  0.0;

    cout << "Enter the number of timesteps " << endl;
    cin  >> N; 

    //
    // compute dt
    //

    dt = (tFinal - tInitial)/double(N);
    cout << "Timestep used " << dt << endl;

    //
    // set initial conditions
    //
    tn = tInitial;
    yn = yInitial;

    //
    // evolve the solution 
    //

    for(i = 1; i <= N; i++)
    {
       yn = EulerStep(yn,tn,dt,f); 
       tn = tn + dt;
    }

    //
    // output results
    //
    cout.precision(14);
    cout << "Solution at t = " << tFinal << " : " << yn << endl;

    return 0;
}

double EulerStep(double yn, double tn, double dt, 
double(*F)(double,double))
{
//
// This routine advances the ordinary differential equation 
//
//                    dy/dt = F(y,t)  
//
// one step using Euler's method. 
//
// The argument double(*F)(double,double) is a function pointer to a double 
// valued function of two double variables that specifies the ordinary differential 
// equation. 
//
    double ynp1;

    ynp1 = yn + dt*F(yn,tn);

    return ynp1;
}