/*
Program RK2.c wrtitten by Erik Cheever
Approximates the solution to a 1st order differential equation using a 2nd order Runge-Kutta

 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 h, float *k1, float *k2);
void main();

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

float xapprox[NDERIVS];				//Intermediate approx of state variables.
int i;								//loop counter.

	k1[0]=3*exp(-4*t)-2*x[0];       //Calculate derivative at beginning of time step.
	
	for(i=0;i<NDERIVS;i++) {		//Approximate state variables at end of time step.
		xapprox[i]=x[i]+h*k1[i];
	}

	k2[0]=3*exp(-4*(t+h))-2*xapprox[0];	//Calculate derivative at end of time step.	
	
}


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

	printf("Welcome to RK2 - 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, h, k1, k2);
		
		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]+k2[i])/2.0;
		}
		
	}
	printf("Program is done.\n");
}