四阶Runge-Kutta法解常微分方程

/**
***四阶Runge-Kutta法***
  
   经典格式：
 y(n+1) = y(n) + h/6 ( K1 + 2*K2 + 2*K3 + K4 )
 K1 = f( x(n) , y(n) )
 K2 = f( x(n+1/2) , y(n) + h/2*K1 )
 K3 = f( x(n+1/2) , y(n) + h/2*K2 )
 K4 = f( x(n+1) , y(n) + h*K3 )

    Runge-Kutta法是基于泰勒展开方法，因而需要所求解具有较好的光滑性。

    属性：差分方法

   《数值分析简明教程》-2 Editon -高等教育出版社- page 105 算法流程图
  
    代码维护：2005.6.14  DragonLord

**/
#include<iostream.h>
#include<stdio.h>
#include<math.h>
/*
举例方程：
 
 y'= y - 2*x / y ( 0<x<1 )
    y(0) = 1

*/
double f(double x,double y)
{
 double re;
 if(x==0)re=1;
 else re=y-2*x/y;
 return re;
}

int main()
{
 double x0,x1,y0,y1,h,k1,k2,k3,k4,y;
 int N;
 while(cin>>x0>>y0>>h>>N)
 {
  int n=0;
 
  for(;n<N;n++)
  {
   x1=x0+h;
   y=sqrt(1+2*x1);
   k1=f(x0,y0);
   k2=f(x0+h/2,y0+h*k1/2);
   k3=f(x0+h/2,y0+h*k2/2);
   k4=f(x1,y0+h*k3);
   y1=y0+h*(k1+2*k2+2*k3+k4)/6;
   
   printf("%.1f %.4f %.4f/n",x1,y1,y);
   x0=x1;
   y0=y1;

  }
 }
 return 0;
}