Решение краевой задачи методом стрельбы

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
 
namespace asdf12345zxc
{
    public partial class Form1 : Form
    {
        double b1 = 0.3;
        public Form1()
        {
            InitializeComponent();
        }
 
        public double f(double x, double y, double z) // y' = z
        {
            return z;
        }
 
        public double g(double x, double y, double z)
        {
            return 2 * Math.Sqrt(x) * z - x * y * (1 - y); // z' = f
        }
 
        public double exSolve(double x) // точное решение
        {
            return -Math.Tan(x);
        }
 
        /**
         * <param name = "a"> Левая граница отрезка </param>
         * <param name="b"> Правая граница отрезка </param>
         * <param name="h"> Шаг по пространству </param>
         * <param name="y0"> Y(a) </param>
         * <param name="yb"> Y(b) </param>
         * <param name="Y"> Массив решений </param>
         * 
         */
        public double[,] rk4(double a, double b, double h, double y0, double yb, double[,] Y)
        {
            double[] k1 = new double[4]; //коэффиценты К
            double[] k2 = new double[4]; //коэффиценты К
            int i = 0;
            Y[0, 0] = y0;
            Y[1, 0] = yb;
            for (double x = a; x < b - h; x += h)
            {
                k1[0] = f(x, Y[0, i], Y[1, i]);
                k2[0] = g(x, Y[0, i], Y[1, i]);
 
                k1[1] = f(x + h / 3, Y[0, i] + k1[0] * (h / 3), Y[1, i] + k2[0] * (h / 3));
                k2[1] = g(x + h / 3, Y[0, i] + k1[0] * (h / 3), Y[1, i] + k2[0] * (h / 3));
 
                k1[2] = f(x + h * 2 / 3, Y[0, i] - k1[1] * (h / 3), Y[1, i] - k2[1] * (h / 3));
                k2[2] = g(x + h * 2 / 3, Y[0, i] - k1[1] * (h / 3), Y[1, i] - k2[1] * (h / 3));
 
                k1[3] = f(x + h, Y[0, i] + k1[0] - k1[1] + k1[2], Y[1, i] + k1[0] - k1[1] + k1[2]);
                k2[3] = g(x + h, Y[0, i] + k1[0] - k1[1] + k1[2], Y[1, i] + k1[0] - k1[1] + k1[2]);
 
                Y[0, i + 1] = Y[0, i] + (h / 8) * (k1[0] + 3 * k1[1] + 3 * k1[2] + k1[3]);
                Y[1, i + 1] = Y[1, i] + (h / 8) * (k2[0] + 3 * k2[1] + 3 * k2[2] + k2[3]);
                i++;
            }
            return Y;
        }
 
        public void toTable(double a, double b, double h, double[,] Y)
        {
            dataGridView1.ColumnCount = 4;
            dataGridView1.RowCount = (int)((b - a) / h + 1);
            dataGridView1.Columns[0].HeaderText = "X";
            dataGridView1.Columns[1].HeaderText = "Метод стрельбы";
            dataGridView1.Columns[3].HeaderText = "Точное решение";
 
            
            dataGridView1_CellContentClick(a,b,h,Y);
        }
        /**
         * <param name = "a"> Левая граница отрезка </param>
         * <param name="b"> Правая граница отрезка </param>
         * <param name="y0"> Y(a) </param>
         * <param name="y1"> Y(b) </param>
         * <param name="eps"> Точность </param>
         * <param name="h"> Размер шага по пространству </param>
         */
        public double[,] shooting_method(double a, double b, double y0, double y1, double h, double eps)
        {
            double nu1 = 1.0, nu2 = 0.8, nu3; // некоторое значение тангенса угла наклона касательной к решению в точке x = a
            int N = (int)((b - a) / h + 1); // количество элементов
            double[,] Y1 = new double[2, N], Y2 = new double[2, N]; // массивы для хранения решений
            Y1 = rk4(a, b, h, y0, nu1, Y1); // решаем систему для nu1
            Y2 = rk4(a, b, h, y0, nu2, Y2); // решаем систему для nu2
            while (Math.Abs(Y2[0, N - 1] - y1) > eps) //
            {
                nu3 = nu2 - ((nu2 - nu1) / (Y2[0, N - 1] - Y1[0, N - 1])) * (Y2[0, N - 1] - y1); // вычисляем новое nu
                nu1 = nu2;
                nu2 = nu3;
                Y1 = Y2;
                Y2 = new double[2, N]; // без этого не работает ?!? о_О
                Y2 = rk4(a, b, h, y0, nu2, Y2); //вычисляем решение для нового nu
            }
            return Y2;
        }
 
        private void Form1_Load(object sender, EventArgs e)
        {
            double a = -1; double b = 1; // отрезок [a; b] 
            double y0 = 0; //y(a)
            double y1 = 1; //y(b)
            double h = 0.09, eps = 0.001; // шаг и точность
            int N = (int)((b - a) / h + 1);
            double[,] Y = shooting_method(a, b, y0, y1, h, eps);
            toTable(a, b, h, Y);
        }
 
        private void dataGridView1_CellContentClick(double a, double b, double h, double[,] Y)
        {
            int i = 0;
            for (double x = a; x <= b; x += h)
            {
                dataGridView1[0, i].Value = x;
                dataGridView1[1, i].Value = Y[0, i];
                dataGridView1[3, i].Value = exSolve(x);
                i++;
            }
        }
    }
}