Смешанная задача для уравнения теплопроводности MATLAB

[OlgaKuz]

Приветствую всех!
Нужно решить смешанную задачу для уравнения теплопроводности в среде MATLAB шеститочечной схемой (вроде как называется "метод Кранка-Николсона"), написала код, который выдаёт фигню, не могу понять, где ошиблась, помогите, пожалуйста!
Мой код:

clear all
close all
clc
%Входные параметры
tau=0.05;
h=0.05;
T=1;
% Сетки
x=0:h:1;
t=0:tau:T;
% Количество точек по X и по T
I=1./h+1;            % по X
J=T/tau+1;           % по T
% Поиск решения
% Находим решение на первом слое (Н.У.)
for i=1:I
    u(1,i)=2.5*tanh(x(i));
end
for j=2:J
    % Прогонка
    alpha(1)=0;
    betta(1)=0;
    alpha(2)=1./(h+1);
    betta(2)=-5.*(tanh(t(j))+1./(cosh(t(j)))^2)/(2.*(h+1));
    for i=2:(I-1)
        % Коэффициенты СЛАУ
        b(i)=-tau/2;
        c(i)=-(h^2+tau);
        a(i)=-tau^2;
        f(i)=-[tau/2.*(u(j-1,i+1)-2.*u(j-1,i)+u(j-1,i-1))+h^2.*u(j-1,i)+5.*(2.*tanh(x(i)-t(j)-0.5.*tau)-1)/(2.*(cosh(x(i)-t(j)-0.5.*tau))^2)]*tau*h^2;
        f(1)=betta(2);
        f(I)=5.*h/(2.*(cosh(t(j)-1))^2);
        a(1)=0;
        c(I)=0;
        % Коэффициенты прогонки
        alpha(i+1)=b(i)/(c(i)-a(i)*alpha(i));
        betta(i+1)=(a(i)*betta(i)+f(i))/(c(i)-a(i)*alpha(i));
    end
    % Поиск решения
    u(j,I)=(f(I)+betta(I))/(1-alpha(I));
    for i=(I-1):-1:1
        u(j,i)=u(j,i+1)*alpha(i+1)+betta(i);
    end
end
% Аналитическое решение
for j=1:J
    for i=1:I
        u0(j,i)=2.5.*tanh(x(i)-t(j));
    end
end
% Построение поверхностей
figure
mesh(x,t,u)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
figure
mesh(x,t,u0)
xlabel('x')
ylabel('t')
zlabel('u0(x,t)')
figure
mesh(x,t,abs(u0-u))
xlabel('x')
ylabel('t')
zlabel('Погрешность')