/*
	Title: Demonstration of trapezoidal runge-kutta method,
		of order two.  It implements the method in formula
		(9.69) of the text, in Section 9.4.
	
	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>
      
double f(double x, double z);
double Y(double x);

int numde;


main()
{
	FILE *fp;
	char q;
	double xzero, yzero, xend, h, x0 , y0;
	double x1, y1, f0, f1, truevalue, error;
	int k, iprint; 
      
	printf("\n Do you want the output on the terminal ?  (y/n)");
	printf("\n If not, then it is stored on file rk2.output : ");
	fflush(stdin);
	q=getchar();
	if ((q == 'N') || (q == 'n'))  
	{
		if ((fp = fopen("rk2.output", "w")) == NULL)  
		{
			printf("\n Cannot open file rk2.output\n");
			return(1);
		}
	}
	else
		fp = stdout;

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

		printf("\n Which Differential equation?");
		printf("\n Give zero to stop : ");	
		scanf("%d", &numde);
		if (numde == 0)  
		{
			fclose(fp);
			return(0);
		}

      	printf("\n Give domain [x0,b] of solution : ");
		scanf("%lf %lf", &xzero, &xend);
		printf("\n What is y0=y(x0)? : ");
		scanf("%lf", &yzero);
		
		while (1)  
		{			
			printf("\n\n Give stepsize h and print parameter iprint");
			printf("\n Let h=0 to try another differential equation : ");
			scanf("%lf %d", &h, &iprint);
			if (h == 0)
				break;

			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  */

			x0 = xzero;	 		/*  Initialize  */
			y0 = yzero;
	
			while (x0 + h <= xend)
			{
				for(k=1; k <= iprint; k++)
				{
					x1 = x0 + h;
					if (x1 > xend)
						break;
					f0 = f(x0, y0);
					f1 = f(x0 + h, y0 + h*f0);
					y1 = y0 + 0.5*h*(f0 + f1);
					x0 = x1;
					y0 = y1;
				}

				/*  Calculate error and print results  */
				if (x1 > xend)
					break;
				truevalue = Y(x1);
				error = truevalue - y1;
				fprintf(fp, "\n x = %10.3e     y(x) = %21.10e", x1, y1);
				fprintf(fp, "     Error = %12.2e", error);
			}
		}
	}
}
      
      
double f(double x, double z)
{
	/*  This defines the right side of
    	the differential equation  */

	const double one = 1.0, two = 2.0;
	double 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);
}


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

	const double one = 1.0, two = 2.0;
	double 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);
}