1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
	
#include<iostream>
#include<stdio.h>
#include<math.h>
 
#define NMAX 500
 
using namespace std;
 
double f1(double x, double y)
{
    // return dy/dx
    return x*y;
}
 
int main()
{
    int i,n;
    double a, b, h;
    double x[NMAX],  y[NMAX];
    double k1, k2, k3, k4;
 
    /* limits */
    cout << "Please enter lower limit, a= " ;
    cin >> a;
    cout << "Please enter lower limit, b= " ;
    cin >> b;
 
    /* discretization */
    cout << "Please enter number of steps: (<=500!) ";
    cin >> n;
    if(n > NMAX)
    {
        cout << "Please try again with fewer steps\n";
        return 0;
    }
    /* initial conditions */
    cout << "Please enter initial value y[0]= ";
    cin >> y[0];
   
    h  = (b - a) / n;      /* step */
    x[0] = a;  /* initial value of x */
 
    /* loop calculations */
    for(i=0; i < n; i++)
    {
        k1 = h * f1(x[i],y[i]);
        k2 = h * f1(x[i] + h/2.0, y[i] + 0.5 * k1);
        k3 = h * f1(x[i] + h/2.0, y[i] + 0.5 * k2);
        k4 = h * f1(x[i] + h, y[i] + k3);
        y[i+1] = y[i] + (k1 + 2.0*k2 +2.0*k3 + k4 ) / 6.0;
        x[i+1] = a + h*(i+1);
    }
 
    /* print resulted table */
    for (i=0; i <= n; i++)
       printf("%5d %10.6f %10.6f\n", i, x[i], y[i]);
   
    return 0;
}