Уравнение нелинейной диффузии Колмогорова-Петрова-Пискунова

% очищаем памяить и экран
clear all
clc 
 
% пределы:
a = 0;
b = 1;
 
k = 0.7; % коэф-т в уравнении
 
syms x t % определяем символьные переменные
u = sin(pi*x)-x; % начальное значение при t=0
%disp('составляем базисные функции и составляем уравнение общего вида')
for i = 1:4
    for j = 1:4
        % базисные функции
        fi(i) = x^(i)*(1-x);
        ro(j) = t^(j)*(1-t);
        C = sym(['C',num2str(i),num2str(j)]); % создаем новую переменную - коэф-т
        u = u + C*fi(i)*ro(j); % общий вид функции, накапливаем сумму
    end
end
 
d2Udx2 = diff(u,'x',2); % вторая производная по х
dUdt = diff(u,'t'); % первая производная по t
 
%disp('Вычисляем невязку')
R = d2Udx2 - dUdt + k*u*(1-u); % невязка
 
II = sym(zeros(1,16)); % пустой массив
%disp('Вычисляем интегралы')
ii = 1;
for i =1:4
    for j = 1:4
        % fprintf('Интеграл №%i ',ii);
        p = expand( R*fi(i)*ro(j) ); % подынтегральная функция (раскрываем скобки)
        % fprintf('Составили подынтегральную ф-цию... ');
        % считаем двойные интегралы
        It = int( collect(p,'t'),'t',a,b); % собираем множители при t и интегрируем по t
        % fprintf('Проинтегрировали по t... ');
        II(ii) = int( collect(It,'x'), 'x', a,b); % собираем множители при х и интегрируем по x
        % fprintf('Проинтегрировали по x \n');
        ii = ii+1; 
    end
end
 
% решаем систему интегралов относительно Сi:
% disp('Решаем систему')
name = symvar(II); % имена коэф-тов
z = solve(II); % решаем систему относительно Cij
 
% подставляем найденный результат в ф-цию общего вида:
% disp('Подставляем в функцию общего вида')
U = subs(u,{name},{z.C11, z.C12, z.C13, z.C14, z.C21, z.C22, z.C23, z.C24, z.C31, z.C32, z.C33, z.C34, z.C41, z.C42, z.C43, z.C44});
 
% рисуем и включаем сетку
disp('Рисуем')
ezsurf(U,[a,b,a,b]), grid on
xlabel('x')
ylabel('t')