/*
implizites und explizites Eulerverfahren für y'(x)=-2000*(y-cos(x))-sin(x)

*/

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

#define ALPHA 200

double f(double x,double y) {
	return -ALPHA*(y-cos(x))-sin(x);
	}
	
double explicit(double x,double y,double h) {
	return y+h*f(x,y);
}
	
double implicit(double x,double y, double h) {
	/* solve yp=y+h*f(xp,yp) for yp
		explicit solution is
		y=yp-h*f(xp,yp)
		 =(1+2000*h)*yp-2000*h*cos(xp)+h*sin(xp)
	*/
	double xp=x+h;
	
	return (y+h*(ALPHA*cos(xp)-sin(xp)))/(1+ALPHA*h);
}


int main(){
	double x,ye,yi,h;
	
	FILE *fp=fopen("imp_euler.dat","w");
	
	x=0; ye=yi=1; h=1.2e-2;
	
	for(;x<1; x+=h){
		fprintf(fp,"%8.5f\t%8.5f\t%8.5f\t%8.5f\n",x,cos(x),ye,yi);
		ye=explicit(x,ye,h);
		yi=implicit(x,yi,h);
	}
	return 0;
}