#include <iostream>
#include <math.h>

using namespace std;

 double f( double x,  double y)
{
        return y/x-y-x;
}
 double norm(const  double *a, const  double *b, int n)
{
        double norma = fabs(b[1] - a[0]);
        for(int i = 0,  j = 1; i < n/2; i++, j+=2)
        {
                norma = norma<fabs(b[j]-a[i]) ? fabs(b[j]-a[i]) : norma;
        }
        return norma;

}

int main()
{
        double e, a,b, *y, *newy, norma,h, y0, k1,k2;
        int n = 10, kst=0;
        cout << "Enter e =";
        cin >> e;
        a = 1;
        b = 4.5;
        y0 = exp(-1) - 1;
        y = new  double [5];
        newy = new  double[10];
        for(int i = 0; i < n; i++) {
                newy[i] = 0;
                y[i/2] = 0;
                }
        do
        {
                delete y;
                y = new  double[n/2];
                for(int i = 0; i < n/2; i++) y[i] = newy[i];
                delete[] newy;
                newy = new  double[n];
                newy[0] = y0;
                h = (b-a)/n;
                for(int i = 1; i < n; i++)
                {
                        k1 = f(a+i*h, newy[i-1]);
                        k2 = f(a+i*h + 2.0/3*h, newy[i-1] + 2.0/3*h*k1);
                        newy[i] = newy[i-1] + h / 4 *k1 + 3.0/4*h*k2;
                }
                norma = norm(y, newy,n);
                n*=2;
                kst++;
        }while(norma>=e);
        n/=2;
        for(int i = 0; i < n; i++)
        {
                cout << a + i*h << " ";
                cout << newy[i] << " ";
                cout << (a + i*h)*(exp(-(a + i*h))-1) << endl;
        }
        cout << b << " ";
        k1 = f(b, newy[n-1]);
        k2 = f(b + 2.0/3*h, newy[n-1] + 2.0/3*h*k1);
        cout << newy[n-1] + h / 4 *k1 + 3.0/4*h*k2 << " ";
        cout << b*(exp(-b)-1) << endl;

        cout << "Iteration = " << kst;

    return 0;
}