//this program implements the regular, modified, and improved Euler's Methods
//	the ODE to be solved is y'=y^2+1 with y(0)=0
//	(see Chapter 5 of DeVries)
//	output will be values of y(x) for the three cases plus exact solution within [0,1]
//	additional outputs are the differences of the numerical solns to the actual soln
//the exact solution of the ODE is y=tan(x)

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

double regEuler(double yy,double xstep);
double modEuler(double yy,double xstep);
double impEuler(double yy,double xstep);
double ode(double yy);

main()
{
	double x;
	double yeuler,ymodified,yimprove,yactual;
	double h;
	
	char filename[25];
	FILE *outfile;

	//front matters
	printf("Euler Methods\n");
	printf("Please enter step size\n");
	scanf("%lg",&h);
	printf("Name for output data file\n");
	scanf("%s",filename);

	outfile=fopen(filename,"w");
	fprintf(outfile,"xx\teuler\tmodeuler\timpeuler\tacteuler\tdiffeu\tdiffmo\tdiffim\n");

	//initialize
	x=0.0;
	yeuler=ymodified=yimprove=yactual=0.0;
	fprintf(outfile,"%20.15lg\t%20.15lg\t%20.15lg\t%20.15lg\t%20.15lg\n",
		x,yeuler,ymodified,yimprove,yactual);
	do
		{
		//increment numerical y
		yeuler=regEuler(yeuler,h);
		ymodified=modEuler(ymodified,h);
		yimprove=impEuler(yimprove,h);
		
		//increment x and actual y
		x+=h;
		yactual=tan(x);

		//output
		fprintf(outfile,"%20.15lg\t%20.15lg\t%20.15lg\t%20.15lg\t%20.15lg\n",
			x,yeuler,ymodified,yimprove,yactual);
		}
	while(x<=1.0);

	return 0;
}

double regEuler(double yy,double xstep)
{
	return yy+xstep*ode(yy);
}

double modEuler(double yy,double xstep)
{
	double ymid;

	ymid=yy+xstep/2.0*ode(yy);
	return yy+xstep*ode(ymid);
}

double impEuler(double yy,double xstep)
{
	double favg;

	favg=(ode(yy)+ode(yy+xstep*ode(yy)))/2.0;
	return yy+xstep*favg;
}

double ode(double yy)//the ode does not depend on x
{
	return yy*yy+1.0;
}