#include<stdio.h>

#include<math.h>

#define f(x,y) (y)/(2.0*(x))  /* dy/dx=f(x,y) */

#define  x0  1.0            /* ¼®·Á         */

#define  y0  1.0            /* y0=y(x0)     */

#define  h   0.1            /* ·»ÞР        */

#define  n   10             /* Y1,¥¥¥,Yn:¶² */

double  x=x0, xi=x0+h, EY=y0, HY=y0, RKY=y0;

void main()

{

    int i;

    double Euler(void);

    double Heun(void);

    double RungeKutta(void);

    printf("\n\n------ Ordinary Differential Equation ------\n\n");

    printf("L.Euler method /  K.Heun method / Runge & Kutta method");

    printf("\n\n\n  dy/dx=f(x,y) , f(x,y)=y/2x \n\n");

    printf(" i    xi  :  Yi(Euler)   Yi(Heun)  Yi(Runge-Kutta)");

    printf("  yi(=sqrt(xi))\n\n");

    printf(" 0  %5.2f :                                       ",x0);

    printf("    %9.6f",y0);

    for(i=1;i<=n;i++)

    {

        printf("\n%2d  %5.2f :",i,xi);

        printf("  %9.6f  %9.6f  %9.6f",Euler(),Heun(),RungeKutta());

        printf("         (%9.6f)",sqrt(xi));

        x=x+h;

        xi=xi+h;

    }

    printf("\n\n");

}

double Euler()

{

    EY=EY+h*f(x,EY);

    return EY;

}

double Heun()

{

    double YY;

    YY=HY+h*f(x,HY);

    HY=HY+0.5*h*(f(x,HY)+f(xi,YY));

    return HY;

}

double RungeKutta()

{

    double A,B,C,D;

    A=h*f(x,RKY);

    B=h*f(x+h*0.5,RKY+A*0.5);

    C=h*f(x+h*0.5,RKY+B*0.5);

    D=h*f(xi,RKY+C);

    RKY=RKY+(A+2.0*B+2.0*C+D)/6.0;

    return RKY;

}