改进Eular方法解常微分方程

/**
***改进Eular方法***

  预报：*Y(n+1)=Y(n) + h * f( x(n) , y(n) );
 
  改进：Y(n+1)=Y(n) + h/2 *［ f( x(n) , y(n) ) + f( x(n+1) , *Y(n+1) ) ］

  h=x(n+1)-x(n)

  平均化形式：
  Yp= Yn + hf( X(n),Y(n) )
  Yc= Yn + hf( X(n+1),Y(p) )
  Y(n+1)= 1/2 * ( Yc+Yp )

  属性：差分方法

 《数值分析简明教程》-2 Editon -高等教育出版社- page 100 算法流程图
  
  代码维护：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,y0,x1,y1,y;
 double yp,yc,h;
 int N;
 while(cin>>x0>>y0>>h>>N)
 {
       
  int n=0;
  for(;n<N;n++)
  {
   x1=x0+h;
   y=sqrt(1+2*x1); //精确值
   
   yp=y0+h*f(x0,y0);
   yc=y0+h*f(x1,yp);
   
   y1=(yp+yc)/2.0;
   printf("%.1f %.4f %.4f/n",x1,y1,y); //输出与精确值比较
   x0=x1;
   y0=y1;
  }
 }
 return 0;
}