/* Draws the graph of approximate solutions of the Cauchy problem for
   a first order ordinary differential equation in the normal form
   by the Euler-Cauchy method.
   Replace the line No.5 by the function in the right hand side
   of the equation to be studied.
   The initial setting is for the solution at x=1 of the equation
   dy/dx = y and for the initial data y0=0 (approximate value of e).
   compile with the following instruction:
   gcc rep8-4.c xgrc.o -lm -lX11 -L/usr/X11R6/lib */
#include<stdio.h>
#include<stdlib.h>
const int IA=0,IB=799,IC=0,ID=599;
double A,B,C,D;
int IGX(double X)
{
    return IA+(X-A)/(B-A)*(IB-IA);
}
int IGY(double Y)
{
    return IC+(Y-C)/(D-C)*(ID-IC);
}
double F(double X,double Y)
{
    return Y;
}
int main(void)
{
    int I,N;
    double H,H2,X,Y,SCH,Y0,YS,U,V,VS;
    char ch;
    init_(IA,IC,IB,ID);
    cls_();
    printf("number of division ? (e.g. 1000) : ");
    scanf("%d",&N);
    printf("range of the variable  X? (e.g. 0, 1) : ");
    scanf("%lf,%lf",&A,&B);
    printf("      of               Y? (e.g. -3, 3) : ");
    scanf("%lf,%lf",&C,&D);
    lcolor_(15);
    line_(IGX(A),IGY(0.0),IGX(B),IGY(0.0));
    line_(IGX(0.0),IGY(C),IGX(0.0),IGY(D));
    lcolor_(9);
    printf("initial value at x=0? (e.g. 1) : ");
    scanf("%lf",&Y);
    H=(B-A)/N;
    H2=H/2;
    printf("A=%lf,B=%lf,Y=%lf",A,B,Y);
    X=A;
    move_(IGX(X),IGY(Y));
    for(I=1;I<=N;I++){
        U=H*F(X,Y);
        YS=Y+U/6;
        U=H*F(X+H2,Y+U/2);
        YS=YS+U/3;
        U=H*F(X+H2,Y+U/2);
        YS=YS+U/3;
        U=H*F(X+H,Y+U);
        Y=YS+U/6;
        X=X+H;
        lineto_(IGX(X),IGY(Y));
    }
    printf("Last value of the solution : %18.15lf",Y);
    fflush(stdout);
    printf("Input q to finish : ");
    do {
       scanf("%c",&ch);
    } while(ch!='q');
    closex_();
    return 0;
}