Здравствуйте. Необходимо построить ВАХ тлеющего разряда на основании системы уравнений, одно из которых является немного изменённым уравнением Роговского.

Параметры, содержащиеся в уравнениях:

Мне в 
теме раздела «Околонаучный софт» написали, что интеграл в правой части второго уравнения не сходится при данных параметрах. Тем не менее все параметры перепроверил, к тому же они соответствуют реально возможным в тлеющем разряде. Есть пособие, в котором приведён график, построенный на основании данных уравнений, но формулы для зависимости 

 не приведено.
Я пытался решить двумя способами:
1. Проинтегрировал-таки правую часть второго уравнения. В итоге одним из слагаемых результата интегрирования был экспоненциальный интеграл. Я решил пренебречь им, поскольку:

Постоянная Эйлера порядка 0,57. Натуральный логарифм 

, где вместо 

 в моём конкретном случае будет 

, мал в силу того, что 

 - протяжённость темного катодного пространства, и выраженное в СИ для моего случая будет составлять от силы 0,5 м. Соответственно сумма ряда тоже не внесёт существенного вклада в подсчёты.
Правильно ли я полагаю, что можно упростить себе задачу и пренебречь экспоненциальным интегралом?
В итоге интеграл в правой части раскрывается, как:

Выразил из первого уравнения 

:

Выделил произведение всех констант, содержащихся в 

 в новую константу:

Левую часть второго уравнения также обозначил за константу, предварительно разделив обе части на 

:

Ввёл новую константу 

:

Поделил обе части на 

. В итоге в правой части осталась только экспонента.
Взял натуральный логарифм от обеих частей уравнения:

Каждую часть обозначил за функцию двух переменных:

Пытался построить графики двух этих функций на одном общем графике поверхности в маткаде. Отдельно маткад их строит, при попытке построить вместе выдаёт ошибку. Но по порядку величин видно, что общих точек у данных функций нет.
Cкрин графика функции 

:

График функции 

:

2. Численный метод. Делал в матлабе. Из первого уравнения выразил 

:

Заменил интеграл на сумму. Подставил 

 во второе уравнение, упростил степень экспоненты:

Вынес из-под знака суммы экспоненту:

Поделил левую часть второго уравнения на эту экспоненту. По сути слева получилась функция от 

. В матлабовском коде ниже она обозначена как 

.
Код из матлаба:
Dif1=1000;% задаём начальное значение, с которым будем сравнивать первый модуль разности
Ukmin=0;%определяем начальное минимальное значение напряжения
i=1;%задаем начальное значение индекса массива данных
for jk=0.000001:0.01:10 %для каждого jk 
    for deltaUk=0.00001:50:3000 %проходимся по всем значениям deltaUk
        %находим dk для определения максимального значения x
        dk=((4)^(1/3)*bi^(1/3)*deltaUk^(2/3)*eps0^(1/3)*(y+1)^(1/3))/(jk^(1/3));
        deltadk=0.01*dk;%определяем шаг изменения x
        Sum1=0;%задаём начальное значение суммы
        %в цикле производим суммирование по x
        for x=0:deltadk:dk
            f1=exp(-((dk^2)/(x-dk)));
            Sum1=Sum1+f1;
        end
        %считаем левую часть уравнения для данного значения deltaUk
        f0=(log(1+(1/y))*exp((B_const*p)/(2*deltaUk)))/(p*A_const);
        %ищем модуль разности правой и левой части
        Dif=abs(f0-Sum1);
        %при решении возникают бесконечности, исключаем их
        if (Dif == inf)
            continue;
        end
        %если разность для следующей итерации оказывается меньше предыдущей
        %заменяем значение Dif1 на новое, равное текущему Dif
        if (Dif<=Dif1)%при выполнении условия
            Dif1=Dif;
            Ukmin=deltaUk;%заменяем минимальное значение напряжения на текущее deltaUk
        end
    end
    %заполняем массив полученными значениями
    Mass(1,i)=jk;
    Mass(2,i)=Ukmin;
    Mass(3,i)=Dif1;
    i=i+1;%увеличиваем значение индекса массива на единицу
end
  Ни один из этих способов не дал адекватных результатов. Разность между правой и левой частью во втором случае составляет 100, при равенстве левой части 0,29. То есть даже близко правая часть на всём диапазоне как 

, так и 

, не приближается по значениям к левой.
Закралось подозрение, что дело в константах, но опять же, как писал выше, перепроверял их, и не один раз.
Если есть кто-то, кто сталкивался с выводом зависимости 

 от 

 для тлеющего разряда, буду рад вашей помощи и советам.
В возможности решения данной задачи не сомневаюсь, поскольку, как уже писал выше, есть график ВАХ, построенный на основе вышеупомянутых уравнений.