/*     Title: Demonstration of adams-bashforth method of order two. 
	   This is from Section 9.5
*/

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

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

double f(double, double);
double Y(double);

int numde;

main()
{
	int kbeg, k, iprint;
	double xzero, xend, yzero, h, truevalue, error;
	double x0, y0, x1, y1, x2, y2, f0, f1;
	

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

		printf("\n Which Differential equation?");
		printf("\n Give zero to stop : ");
		scanf("%d", &numde);
		if (numde == 0)
			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;

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

			/*  Initialize  */
			x0 = xzero;	
			y0 = yzero;
			f0 = f(x0,y0);
			x1 = x0 + h;
			y1 = y0 + h*f0;
			f1 = f(x1,y1);
			if (iprint > 1) 
				kbeg = 2;
			else  
			{
				kbeg = 1;
				x2 = x1;
				y2 = y1;
				truevalue = Y(x2);
				error = truevalue - y2;
				printf("\n x = %10.3e     y(x) = %17.10e", x2, y2);
				printf("     Error = %9.2e", error);
			}
			
			/*  Begin main loop for computing solution of 
				differential equation  */
	
			while (x0 + two*h <= xend) 
			{
				for(k=kbeg; k <= iprint; k++) 
				{
					kbeg = 1;
					x2 = x1 + h;
					if (x2 > xend)
						break;
					y2 = y1 + h*(three*f1 - f0)/two;
					x0 = x1;
					x1 = x2;
					y0 = y1;
					y1 = y2;
					f0 = f1;
					f1 = f(x1,y1);
				}
			
				/*  Calculate error and print results  */

				if (x2 > xend)
					break;
				truevalue = Y(x2);
				error = truevalue - y2;
				printf("\n x = %10.3e     y(x) = %19.10e", x2, y2);
				printf("    Error = %11.2e", error);
			}
		}
	}
}
      

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

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

	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);
}