/*
/*編程者：劉雯
/*使用經典四階龍格-庫塔算法進行高精度求解　
/* 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 fu(float x1,float x2, float c)//控制函數u,k取50，e取1
{    float s,u1,u2;
       s=x1+c*x2;
	   if(s>0)
		   {u1=(1+50*c*x1+47*x2+c*x2)/9;
		    return u1;}
	   else 
	   {if(s<0){u2=(-1+50*c*x1+47*x2+c*x2)/9;
		      return u2;}
	   else {printf("Error"); return 0;}
	  
	   }


}

double f2(float mid_value,float u)  //狀態變量x2的方程
{
	double rus;//,pi=3.1415926
	rus=-3*mid_value-9*u;//-11.25*sin(pi*i/20)+22.5*cos(pi*i/20)
	return(rus);
}
float f1(float mid_value)//狀態變量x1的方程
{
    float rus;
    rus=mid_value;
    return(rus);
}
void main()
{
	float h,Value_u,c;  //需要的變量 
	int n;    //計算出的點的個數
 	float k1,k2,k3,k4,k5,k6,k7,k8;
    float x1[40],x2[40]; //用于存放計算出的常微分方程數值解
	int i=0;
	

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

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

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

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

// 以下為核心程序
    //n=((x1-x0)/h);
	n=40;
	
    for(i=0;i<n;i++)
	{  
	  Value_u=fu(x1[i],x2[i],c);
//	  Value_x2=f2(x2[i],Value_u);
//      Value_x1=f1(x2[i]);
      printf("Value_u=%f\n",Value_u);
      //求x2
      k1=f2(x2[i],Value_u); //求K1
      k2=f2((x2[i]+k1*h/2),Value_u); //求K2
      k3=f2((x2[i]+k2*h/2),Value_u); //求K3
      k4=f2((x2[i]+k3*h),Value_u); //求K4
 
	  x2[i+1]=x2[i]+(k1+2*k2+2*k3+k4)*h/6; //求x2[i+1]
  //    printf("k1=%f\n",k1);
 //     printf("k2=%f\n",k2);
 //     printf("k3=%f\n",k3);
 //     printf("k4=%f\n",k4);
      printf("x2[%i]=%f\n",i+1,x2[i+1]); 

      //求x1
      k5=f1(x2[i]); //求K5
      k6=f1(x2[i]+k5*h/2); //求K6
      k7=f1(x2[i]+k6*h/2); //求K7
      k8=f1(x2[i]+k7*h); //求K8
 
	  x1[i+1]=x1[i]+(k5+2*k6+2*k7+k8)*h/6; //求x1[i+1]
      
//	  printf("k5=%f\n",k5);
//      printf("k6=%f\n",k6);
//      printf("k7=%f\n",k7);
//      printf("k8=%f\n",k8);
      printf("x1[%i]=%f\n",i+1,x1[i+1]);
      
	 
	}
 
}