close all
clear all
clc

%% data
l = 3;
T = 3;
nx = 100;
nt = 150;
hx = l/(nx-1);
ht = T/nt;


%% equation's data
k = @(u) u.^2;
m = @(t) 1 + t^0.5;
c = @(u) 1 + 0.1 * u;
r = @(u) 1 + 0.05 * u.^3;
f = @(u) 2 * u;


Y0 = zeros(nx, 1);
RES = zeros(nx, nt);
RES(:,1) = Y0;
%% solver
eps = 1e-3;
ax = axes;
count = 1;
x = linspace(0, l, nx);
for j = 2:nt
    Yt = Y0;
    while 1
        Ys = Y0;
        K = k(Ys);
        C = c(Ys);
        R = r(Ys);
        F = f(Ys);
        Y0 = makeData(hx, ht, nx, F, K, C, R, m(ht * (j - 1)), Yt);
        if max(abs(Y0 - Ys)) / max(abs(Ys)) <= eps
            break;
        end
    end
    RES(:, j) = Y0;
    if mod(j, 0.2*T/ht + 1) == 0
        axi = subplot(2,2,count)
        plot(axi,x,Y0,'b')
        title(['T = ', num2str(ht*(j - 1))])
%         plot(x, Y0)
        xlim(axi, [0 (T + 0.5)]);
        ylim(axi, [0 (T + 0.5)]);
        drawnow;
        count = count + 1;
    end
end