/*
/*編程者：劉艮平
/*完成日期：2004.5.29
/*e.g:y'=y-2x/y,x∈[0,0.6]
/*　　y(0)=1
/*使用經典四階龍格-庫塔算法進行高精度求解　
/* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6
/* K1=h*f(xi,yi)
/* K2=h*f(xi+h/2,yi+K1/2)
/* K3=h*f(xi+h/2,yi+K2/2)
/* K4=h*f(xi+h,yi+K3)
*/
#include <stdio.h>
#include <math.h>

float f(float den,float p0)  //要求解的微分方程的右部的函數 e.g: y-2x/y
{
	float rus;
 //   den=w/(W0+sl);
 //   rus=k*A*f/Wp*sqrt(RTp)-(k-1)*.....
	rus=p0-2*den/p0;
	return(rus);
}

//float fden()
//{

//}


void main()
{
	float x0; //范圍上限
	int x1;   //范圍下限
	float h;  //步長
	int n;    //計算出的點的個數
 	float k1,k2,k3,k4;
    float y[3]; //用于存放計算出的常微分方程數值解
	int i=0;
	int j;

//以下為函數的接口
	printf("intput x0:");
    scanf("%f",&x0);

	printf("input x1:");
	scanf("%f",&x1);

    printf("input y[0]:");
    scanf("%f",&y[0]);

	printf("input h:");
    scanf("%f",&h);

// 以下為核心程序
    n=((x1-x0)/h);
	n=3;
	
    for(j=0;j<n;j++)
	{   

      k1=h*f(x0,y[i]); //求K1
      k2=h*f((x0+h/2),(y[i]+k1/2)); //求K2
      k3=h*f((x0+h/2),(y[i]+k2/2)); //求K3
      k4=h*f((x0+h),(y[i]+k3)); //求K4
 
	  y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求yi+1
	  x0+=float(0.2);
      printf("y[%f]=%f\n",x0,y[i+1]); 
      ++i;
	 
	}
 
}