/*
      Title: Demonstration of fehl45.
 
      This solves the initial value problem
         y'(x) = f(x,y(x)) ,  x0 <= x <= b ,  y(x0)=y0.
      The function f(x,z) is defined below, along with the true
      solution y(x).  The number of the problem to be solved
      is specified by the input variable 'numde', which is used
      in the functions 'f' and 'y'.  The program will request
      the problem parameters, along with the values of 'h' and
      'iprint'.  'h' is the stepsize, and 'iprint' is the number
      of steps between each printing of the solution.
      Use h=0 and numde=0 to stop the program.
*/


#include <stdio.h>
#include <math.h>

const float zero = 0.0, one = 1.0, two = 2.0, three = 3.0;

float f(float x, float z);
float Y(float x);
void fehl45(float (*f)(float, float), float *y, float *x, 
		float h, float *trun);

int numde;

main()
{
	int k, iprint;
	char q;
	float xzero, xend, yzero, h, truevalue, error;
	float x0, y0, x1, trun;
	FILE *fp;

	
	printf("\n Do you want the output on the terminal ?  (y/n)");
	printf("\n If not, then it is stored on file rkf45.output : ");
	fflush(stdin);
	q=getchar();

	if ((q == 'N') || (q == 'n'))  
	{
		if ((fp = fopen("rkf45.output", "w")) == NULL)  
		{
			printf("\n Cannot open file rkf45.output\n");
			return(1);
		}
	}
	else
		fp = stdout;
	

	while (1)	
	{
		
		/* Input problem parameters */

		printf("\nWhich Differential equation?");
		printf("\nGive zero to stop : ");	
		scanf("%d", &numde);
		if (numde == 0) 
		{
			fclose(fp);
			return(0);
		}
      		
		printf("\nGive domain [x0,b] of solution : ");
		scanf("%f %f", &xzero, &xend);
		printf("\nWhat is y0=y(x0)? : ");
		scanf("%f", &yzero);

		while (1)  
		{
			printf("\n\nGive stepsize h and print parameter iprint");
			printf("\nLet h=0 to try another differential equation : ");
			scanf("%f %d", &h, &iprint);
			if (h == 0)
				break;

			/*  Initialize  */

			x0 = xzero;
			y0 = yzero;
      		
			fprintf(fp, "\n\n Equation%2d     xzero = %9.2e", numde, xzero);
			fprintf(fp, "     b = %9.2e     yzero = %12.5e", xend, yzero);
			fprintf(fp, "\n Stepsize = %10.3e     ", h);
			fprintf(fp, "Print Parameter = %3d\n", iprint);
      		

			/*  Begin main loop for computing solution 
				of differential equation  */
			
			while (x0 + h <= xend)
			{
				
				for (k=1; k <= iprint; k++)  
				{
					x1 = x0 + h;
					if (x1 > xend) 
						break;
					fehl45(f,&y0,&x0,h,&trun);
				}
					/*  Calculate error and print results  */

				if (x1 > xend) 
					break;
				truevalue = Y(x0);
				error = truevalue - y0;
				fprintf(fp, "\nx = %8.2e   y(x) = %15.6e", x0, y0);
				fprintf(fp, "  Error = %11.2e  trun = %11.2e", error, trun);
			}
		}
	}
}


void fehl45(float (*f)(float, float), float *y, float *x, float h, float *trun)
{ 
/* 
	This computes a single step of the runge-kutta-fehlberg method
	of order (4,5) for solving the differential equation
       	y'(x)=f(x,y(x)).

	Input:
	(1) 'f' is the name of a function subprogram for evaluating the
    	derivative f(x,y) that defines the differential equation.
	(2) y is the known value of the solution at the point x.
	(3) h is the stepsize in solving the differential equation.
    	The solution is to be found at x+h.
	
	Output:
	(1) x will be replaced by x+h.
	(2) y will be  replaced by the numerical solution based on
    	using the fourth order member of the fehlberg formula
    	at x+h.
	(3) trun will be set to the estimated truncation error in the
    	fourth order solution y; it is obtained using the fifth
    	order fehlberg solution.
	
	To solve the differential equation on a longer interval,
	call this routine repeatedly.
*/
	

      	
	/* The following constants are needed in evaluating
      	the fehlberg formula  */

	const float b30 = 1932.0/2197.0,b31=-7200.0/2197.0,
				b32=7296.0/2197.0,b40=439.0/216.0,b42=3680.0/513.0,
     			b43=-845.0/4104.0, b50 = -8.0/27.0, b52=-3544.0/2565.0,
     			b53=1859.0/4104.0,b54=-11.0/40.0,a3=12.0/13.0,
     			c0=25.0/216.0,c2=1408.0/2565.0,c3=2197.0/4104.0,
     			d0=16.0/135.0,d2=6656.0/12825.0,d3=28561.0/56430.0,
     			d5=2.0/55.0;

	/*  The following arrays are the constants in the fehlberg formulas  */
      	
	float a[6] = {0.0,0.25,0.375,a3,1.0, 0.5};
	float b[6][5] = {	{0.0, 0.0, 0.0, 0.0, 0.0},
						{0.25, 0.0, 0.0, 0.0, 0.0},
			  			{0.09375, 0.28125, 0.0, 0.0, 0.0},
			  			{b30, b31, b32, 0.0, 0.0},
			  			{b40, -8.0, b42, b43, 0.0},
			  			{b50, 2.0, b52, b53, b54} };
	float c[5] = {c0, 0.0, c2, c3, -0.2};
	float d[6] = {d0, 0.0, d2, d3, -0.18, d5};
	float v[6]; 
	int i,j;
	float fsum, sum, y4, y5;
	

	/*  Evaluate the derivatives  */

	v[0] = (*f)(*x,*y); 	
	for (i=1; i <= 5; i++)	
	{
		sum  = zero;
		for (j=0; j <= i-1; j++) 
			sum = sum + b[i][j]*v[j];
		v[i] = (*f)((*x) + a[i]*h, (*y) + h*sum);
	}	

	/* Evaluate the fourth and fifth order formulas */

	fsum = zero;
	for(i=0; i <= 4; i++)
		fsum = fsum + c[i]*v[i];
	y4 = *y + h*fsum;
	fsum = zero;
	for(i=0; i <= 5; i++)
		fsum = fsum + d[i]*v[i];
	y5 = *y + h*fsum;

	/*  Set the output variables  */

	*x = *x + h; 		
	*y = y4;
	*trun = y5 - y4;
}
	

float f(float x, float z)
{
	/*  This defines the right side of
    	the differential equation  */

	float result;

	switch (numde)  
	{
	case 1: result = -z;
			break;
	case 2: result = (z + x*x - two)/(x + one);
			break;
	case 3: result = -z + two*cos(x);
			break;
	}
	
	return(result);
}


float Y(float x)
{
	/*  This gives the true solution of
    	the initial value problem. */

	float result, x1;

	switch (numde)  
	{
	case 1: result = exp(-x);
			break;
	case 2: x1 = x + one;
			result = x*x - two*(x1*log(x1) - x1);
			break;
	case 3: result = sin(x)  +cos(x);
			break;
	}

	return(result);
}