/*
	Title: Demo of trapezoidal runge-kutta method for systems.
		It implements (9.139)-(9.140) from the text.

	This solves the initial value problem for the first order
	system
   	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 subroutines 'fcn' and 'TrueSoln'.  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.
	
	To increase the order of the systems which can be
	handled, increase the size of 'MAXN' in the parameter
	statement given below.
*/

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

#define	MAXN	3

void fcn(double x,double z[MAXN+1],double f[MAXN+1]); 
void TrueSoln(double x,double y[MAXN+1]);

int numde;


main()
{
	const double zero = 0.0;
	
	double y0[MAXN + 1], y1[MAXN + 1], y[MAXN + 1], yzero[MAXN + 1];
	double f0[MAXN + 1], f1[MAXN + 1], h, x0, x1, xend, xzero; 
	int j, k, iprint, n;
	FILE *fp;


	if ((fp = fopen("rk2.output", "w")) == NULL)  
	{
		printf("\n Cannot open file rk2.output\n");
		return(1);
	}
	
	while (1)  
	{
		
		/*  Input problem parameters  */
		
		printf("\n\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 Give the order of the system.  It should be");
		printf(" less than or equal to MAXN = %d : ", MAXN);
		scanf("%d", &n);

		printf("\n Give the components of y0, one at a time ");
		for (j=1; j <= n; j++)  
		{
			printf("\n Give component : ");
			scanf("%lf", &yzero[j]);
		}
				
		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;
							
			/*  Initialize  */
			x0 = xzero;
			for (j=1; j <= n; j++)
				y0[j] = yzero[j]; 

			fprintf(fp, "\n\n\n Equation %d     xzero = %9.2e", numde, xzero);
			fprintf(fp, "     b = %9.2e\n Stepsize = %10.3e", xend, h);
			fprintf(fp, "     Print Parameter = %3d", iprint);
			fprintf(fp, "\n Order = %1d\n\n j    y0[j]", n);
			for(j=1; j <= n; j++)
				fprintf(fp, "\n%2d  %9.2e", j, yzero[j]);
			printf("\n\n          x     j        approx y[j]          Error");
			
			/*  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;
					fcn(x0,y0,f0);
					for (j=1; j <= n; j++)  
						y1[j] = y0[j] + h*f0[j];
					fcn(x1,y1,f1);
					for (j=1; j <= n; j++)  {
						y1[j] = y0[j] + 0.5*h*(f0[j] + f1[j]);
						y0[j] = y1[j];
					}
					x0 = x1;
				}
						
				/*  Calculate error and print results  */
				if (x1 > xend)
					break;
				TrueSoln(x1,y);
				for (j=1; j <= n; j++)  
				{
					printf("\n %12.5f   %1d     %15.6e", x1, j, y1[j]);
					printf("    %12.2e", y[j] - y1[j]);
				}
			}
		}
	}
}
 
      
void fcn(double x,double z[MAXN+1],double f[MAXN+1]) 
{
	/* This defines the right side of the
	    differential equation, f(x,z)  */
      
	const double two = 2.0, three = 3.0, four = 4.0, five = 5.0;
	double cs, sn;
      
	switch (numde)  
	{
	case 1: cs = cos(x);
      		sn = sin(x);
      		f[1]= z[1] - two*z[2] + four*cs - two*sn;
      		f[2] = three*z[1] - four*z[2] + five*(cs - sn);
      		break;
	case 2: f[1] = z[2];
      		f[2] = z[3];
      		f[3] = -three*(z[3] + z[2]) - z[1] - four*sin(x);
      		break;
	}
}

     
 
void TrueSoln(double x,double y[MAXN+1])
{

	/* This gives the true solution of
	   the initial value problem, y(x)  */
      
	const double two = 2.0;
	double sn, cs;
      
	switch (numde)  
	{
	case 1: cs = cos(x);
      		sn = sin(x);
      		y[1] = sn + cs;
      		y[2] = two*cs;
      		break;
	case 2: sn = sin(x);
      		cs = cos(x);
      		y[1] = sn + cs;
      		y[2] = cs - sn;
      		y[3] = -sn - cs;
      		break;
	}
}