#include<iostream.h>
#include<math.h>
#include<iomanip.h>

float h;
double f0(float );
double f(double x,double y);
double Euler(double x,double y);
double ModEuler(double x,double y);
double RK(double x,double y);

void main(){
	int X=2,n;
    double *x,*y,*y1,*y2,*y3;
	cout<<"請選擇步長:0.1    0.2    0.4"<<endl;
	cout<<"h=";
	cin>>h;
	n=X/h+2;
    x=new double[n];
	y=new double[n];
	y1=new double[n];
	y2=new double[n];
	y3=new double[n];
	x[0]=0;
	for(int i=1;i<n;i++)x[i]=x[i-1]+h;	
	for( i=0;i<n;i++)y[i]=f0(x[i]);
	y1[0]=y[0],y2[0]=y[0],y3[0]=y[0] ;
	for( i=1;i<n;i++)y1[i]=Euler(x[i-1],y1[i-1]);
	for( i=1;i<n;i++)y2[i]=ModEuler(x[i-1],y2[i-1]);
	for( i=1;i<n;i++)y3[i]=RK(x[i-1],y3[i-1]);
	cout<<"x        "<<"精確解   "<<"歐拉方法  "<<"改進的歐拉方法"<<" 經(jīng)典RK方法"<<endl;
   	for( i=0;i<n;i++)
		cout<<setw(3)<<x[i]<<setw(12)<<y[i]<<setw(11)<<y1[i]<<setw(15)<<y2[i]<<setw(12)<<y3[i]<<endl;
    delete[]x;
	delete[]y;
	delete[]y1;
	delete[]y2;
	delete[]y3;
}

double f0(float x){
	return sqrt(4+5*exp(-x*x));
}
double f(double x,double y){
	return 4*x/y-x*y;
}
double Euler(double x,double y){
	return y+h*f(x,y);
}
double ModEuler(double x,double y){
	double p,c;
	p=y+h*f(x,y);
	c=y+h*f(x+h,p);
	return 0.5*(p+c);
}
double RK(double x,double y){
	double k1,k2,k3,k4;
	k1=f(x,y);
	k2=f(x+h/2,y+h/2*k1);	
	k3=f(x+h/2,y+h/2*k2);		
	k4=f(x+h,y+h*k3);
	return y+h/6*(k1+2*k2+2*k3+k4);	
}