using System;

namespace PredictorCorrectorSolver
{
    class PredictorCorrector
    {
        public static double[] x;
        public static double[] y;
        public static double[] f;

        public static void PredCorrSolve(double x0, double xn, double y0, double h)
        {
            int n = (int)(((xn - x0) / h) + 1.0);
            x = new double[n]; x[0] = x0;
            y = new double[n]; y[0] = y0;
            f = new double[n];
            for (int i = 1; i < n; i++)
            {
                x[i] = x[0] + i * h;
            }

            for (int i = 0; i <= 2; i++)
            {
                y[i + 1] = RungeKutta(x[i], y[i], h);
                f[i] = Function(x[i], y[i]);
            }

            for (int i = 3; i < (n - 1); i++)
            {
                f[i] = Function(x[i], y[i]);
                y[i + 1] = y[i] + (h / 24) * (55 * f[i] - 59 *
                    f[i - 1] + 37 * f[i - 2] - 9 * f[i - 3]);
                double tmp;
                while (true)
                {
                    f[i + 1] = Function(x[i + 1], y[i + 1]);
                    tmp = y[i] + (h / 24) * (9 * f[i + 1] + 19 *
                        f[i] - 5 * f[i - 1] + f[i - 2]);
                    if (!(Math.Abs(tmp - y[i + 1]) < 0.00001))
                        y[i + 1] = tmp;
                    else break;
                }
            }
        }

        private static double RungeKutta(double xi, double yi, double h)
        {
            double[] k = new double[4];
            k[0] = Function(xi, yi);
            k[1] = Function(xi + h / 2, yi + h * k[0] / 2);
            k[2] = Function(xi + h / 2, yi + h * k[1] / 2);
            k[3] = Function(xi + h, yi + h * k[2]);
            return yi + (h / 6) * (k[0] + 2 * k[1] + 2 * k[2] + k[3]);
        }

        private static double Function(double xi, double yi)
        {
            return 2 * (Math.Pow(xi, 2) + yi);
        }
    }
}