/*
Program Euler1.c wrtitten by Erik Cheever.
Approximates the solution to a 1st order differential equation using Eulers's method

 The problem to be solved is:

y'(t)+2*y(t)=3*exp(-4*t)

Note: this problem has a known exact solution
 y(t)=2.5*exp(-2*t)-1.5*exp(-4*t)
 
This program is quite general, it is easily adapted to almost any problem
by changing the routine "CalcDerivs" that calculates the derivatives, and the value
of NDERIVS.
*/

#include <stdio.h>		//I/O library
#include <math.h>		//Math library

#define	NDERIVS	1		//Number of derivatives


void CalcDerivs(int numDerivs, float *x, float t, float *k1);
void main();

void CalcDerivs(int numDerivs, float *x, float t, float *k1) {
// numDerivs is the number of derivatives to be calculated
// x is an array containing the state variables
// t is time
// k1 is the array of calculated derivatives.

	k1[0]=3*exp(-4*t)-2*x[0];        //Calculate derivative;
}


void main() {
	float yexact, ystar;			//Exact and approximate values of y
	float x[NDERIVS], k1[NDERIVS];	//Arrays to hold state variables and derivatives.
	float h=0.1, t=0.0, tmax=4.0;	//Step size, time, final time.
	int i;							//loop counter
	

	printf("Welcome to Euler - the ODE solver\n");
	
	x[0]=1;							//Initial conditions.
	ystar=x[0];						//output is our only state variable in this case.
	
	
	for (t=0;t<tmax;t=t+h) {		//for loop to iterate through time.
		CalcDerivs(NDERIVS, x, t, k1);
		
		ystar=x[0];					//Get approx and exact y.
		yexact=2.5*exp(-2*t)-1.5*exp(-4*t);	

		//Print results.
		printf("t=%4.2f, yexact=%6.4f, ystar=%6.4f\n",t,yexact,ystar);
		
		for(i=0;i<NDERIVS;i++) {	//Update state variables.
			x[i]=x[i]+h*k1[i];
		}
		

	}

	printf("Program is done.\n");
}