Delay DE Max

        double[,] alpha;
        double[,] beta;
        int k_max;

        int T0 = 3, T = 15, h = 3;  // k in T0, T

        // max. element from alpha(beta)[n,k1] to alpha(beta)[n,k2]
    public double MAX_S(double [,] arr, int n, int k1, int k2)
        {
            double pmax = arr[n, k1];

            k_max = k1;  // position of the maximal element
            for (int s = k1+1; s <= k2; s++)
                if (pmax < arr[n, s])
        {
            pmax = arr[n, s]; k_max = s;
        }
            return pmax;
        }

        private void CalcAlphaBeta()
        {
            alpha = new double[10, 20];
            beta  = new double[10, 20];

            int N, k;
            double x, max_a, max_b;

            for (k = T0 - h; k <= T; k++)
            {
                alpha[0, k] = -k;
                beta[0, k]  = k;
            }

            double Ln, Mn, Kn, Pn;
            double na, nb;
            double q1, q4, q3, q2, q5, q6, Q1, Q4, Q3, Q2, Q5, Q6;
            double max_an, max_bn, Pa, Pb;
            int R;
            double n2 = 1;

            for (N = 1; N < 9; N++)
            {
                Ln = Math.Min(-alpha[N - 1, 0], 0.50 - alpha[N - 1, 1]); Kn = Math.Min(Ln, n2);
                alpha[N, 0] = -Kn * Ln;
                alpha[N, 1] = 0.5 - Kn * Ln;
                alpha[N, 2] = Math.Pow(1.1, 0.1 * alpha[N - 1, 15]);

                Mn = Math.Min(beta[N - 1, 0], beta[N - 1, 1] - 0.50); Pn = Math.Min(Mn, n2);
                beta[N, 0] = Pn * Mn;
                beta[N, 1] = 0.5 + Pn * Mn;
                beta[N, 2] = Math.Pow(1.1, 0.1 * beta[N - 1, 15]);

                na = 4;
                for (k = T0; k <= T; k++)
                {
                    max_a  = MAX_S(alpha, N-1, k-h, k);
                    max_b  = MAX_S(beta, N-1, k-h, k);
                    max_an = MAX_S(alpha, N, k-h, k-1);
                    max_bn = MAX_S(beta, N, k-h, k-1);

                    if (k == 3)
                        R = k - 1;
                    else if (k <= 8)
                        R = k - 2;
                    else
                        R = k - 3;

                    nb = na * 4;

                    q1 = nb - alpha[N - 1, k];
                    q3 = nb - max_a;
                    q5 = nb - alpha[N - 1, R];

                    q2 = nb - beta[N - 1, k];
                    q4 = nb - max_b;
                    q6 = nb - beta[N - 1, R];

                    Q1 = na / q1 / q1;  Q3 = na / q3 / q3;  Q5 = na / q5 / q5;
                    Q2 = na / q2 / q2;  Q4 = na / q4 / q4;  Q6 = na / q6 / q6;

                    Pa = na / q1 + na / 2 / q3 + na / q5 - alpha[N - 1, k] * Q1 - max_a * Q3 - alpha[N - 1, R] * Q5;
                    x = (Pa + alpha[N, k - 1] - alpha[N, R] * Q5) / (1 - Q1 - Q3);
                    if (x > max_an)
                        alpha[N, k] = x;
                    else
                        alpha[N, k] = (Pa + alpha[N, k - 1] - alpha[N, R] * Q5 + Q3 * max_an) / (1 - Q1);

                    Pb = na / q2 + na / 2 / q4 + na / q6 - beta[N - 1, k] * Q2 - max_b * Q4 - beta[N - 1, R] * Q6;
                    x = (Pb + beta[N, k - 1] - beta[N, R] * Q6) / (1 - Q2 - Q4);
                    if (x > max_bn)
                        beta[N, k] = x;
                    else
                        beta[N, k] = (Pb + beta[N, k - 1] - beta[N, R] * Q6 + Q4 * max_bn) / (1 - Q2);

                    na = na * 2;
                }
                n2 = n2 / 2;
            }
        // the result is in the arrays 'alpha' and 'beta'
        }