clear all
% This code is written by Chiu-Yen Kao for Math 865L Spring,2009
% for tumor virotherapy model by Jin Wang and Jianjun Tian


N = 41; % number of point in x
dx = 1/(N-1); 
xx = (0:(N-1))'*dx;
dt = 0.5*(dx)^2
max_iter = 10001;
t = (0:(max_iter-1))*dt;
R = zeros(1,max_iter);
U = zeros(N,1);X = zeros(N,1);Y = zeros(N,1);Z = zeros(N,1);V = zeros(N,1);
U1 = zeros(N,1);X1 = zeros(N,1);Y1 = zeros(N,1);Z1 = zeros(N,1);V1 = zeros(N,1);
U2 = zeros(N,1);X2 = zeros(N,1);Y2 = zeros(N,1);Z2 = zeros(N,1);V2 = zeros(N,1);
lambda = 2.0;
burst = 130;
beta = 0.07*burst;
k = 2.0;
s = 56;
w = 20;
D = 3.6;
delta = 5.6;
k0 = 1.0;
gamma = 2.5;
mu = 2.1;

% initialization
R(1) = 2; % mm at t=0


% first iteration
iter = 1;
X(:) = 0.84;
Y(:) = 0.10;
Z(:) = 0.06;
a = 2.1812;
V(:) = (a*exp(-4*xx.^2/a^2));
test = sum(V(:).*xx.^2*dx)-0.5*V(N)*dx;

U(1) = 0;
for i = 1:N-1
        Fp = lambda*X(i+1)-mu*(1-X(i+1)-Y(i+1)-Z(i+1))+...
            s*Y(i+1)*Z(i+1)-w*Z(i+1)^2;
        Fm = lambda*X(i)-mu*(1-X(i)-Y(i)-Z(i))+...
            s*Y(i)*Z(i)-w*Z(i)^2;
        U(i+1) = (xx(i)^2*U(i)+R(iter)/2*dx*(xx(i+1)^2*Fp+xx(i)^2*Fm))/(xx(i+1)^2);
end

R(iter+1) = R(iter)+dt*(U(N));

for i = 1:N
    A = (U(i)-xx(i)*U(N))/R(iter+1);
    F = lambda*X(i)-mu*(1-X(i)-Y(i)-Z(i))+s*Y(i)*Z(i)-w*Z(i)^2;
    if i==1 | i==N
        X1(i) = X(i)+ dt*(lambda*X(i)-beta*V(i)*X(i)-F*X(i));
        Y1(i) = Y(i)+ dt*(beta*V(i)*X(i)-k*Y(i)*Z(i)-delta*Y(i)-F*Y(i));
        Z1(i) = Z(i)+ dt*(s*Y(i)*Z(i)-w*Z(i)*Z(i)-F*Z(i));
    else
        X1(i) = X(i) - dt*A*(X(i)-X(i-1))/dx+...
        dt*(lambda*X(i)-beta*V(i)*X(i)-F*X(i));
        Y1(i) = Y(i) - dt*A*(Y(i)-Y(i-1))/dx+...
        dt*(beta*V(i)*X(i)-k*Y(i)*Z(i)-delta*Y(i)-F*Y(i));
        Z1(i) = Z(i) - dt*A*(Z(i)-Z(i-1))/dx+...
        dt*(s*Y(i)*Z(i)-w*Z(i)*Z(i)-F*Z(i));
    end
end

U1(1) = 0;
for i = 1:N-1
        Fp = lambda*X1(i+1)-mu*(1-X1(i+1)-Y1(i+1)-Z1(i+1))+...
            s*Y1(i+1)*Z1(i+1)-w*Z1(i+1)^2;
        Fm = lambda*X1(i)-mu*(1-X1(i)-Y1(i)-Z1(i))+...
            s*Y1(i)*Z1(i)-w*Z1(i)^2;
        U1(i+1) = (xx(i)^2*U1(i)+R(iter+1)/2*dx*(xx(i+1)^2*Fp+xx(i)^2*Fm))/(xx(i+1)^2);
end


A1 = -(xx.*U1(N)/R(iter+1)+2*D./(R(iter+1)^2)./xx);
A1(1) = 0;
A2 = -D./R(iter+1)^2;
S = 2*V+2*delta*dt*Y1;
L1 = 2*dt/dx/dx*(A2)-dt/dx*A1;
D1 = 2-4*dt/dx^2*A2+2*k0*dt*Z1+2*gamma*dt;
UU1 = 2*dt/dx^2*A2+dt/dx*A1;
L1(N) = 4*dt/dx/dx*(A2)-dt/dx*A1(N);
UU1(1) = 4*dt/dx^2*A2+dt/dx*A1(1);
V1 = tridiagSolve(L1(2:N), D1(1:N), UU1(1:N-1), S);


for iter = 2:max_iter-1
    if mod(iter,20)==0
    plot(xx,X1,xx,Y1,xx,Z1,xx,U1,xx,V1)
    legend('X','Y','Z','U','V')
    title(['T = ' num2str((iter-1)*dt) '    R = ' num2str(R(iter))]);
     drawnow
    end
    
    
    R(iter+1) = R(iter)+0.5*dt*(3*U1(N)-U(N));
    for i = 1:N
        A = (U1(i)-xx(i)*U1(N))/R(iter+1);
        F = lambda*X1(i)-mu*(1-X1(i)-Y1(i)-Z1(i))+s*Y1(i)*Z1(i)-w*Z1(i)^2;
        if (i==1) | (i==N)
        X2(i) = X(i)+ 2*dt*(lambda*X1(i)-beta*V1(i)*X1(i)-F*X1(i));
        Y2(i) = Y(i)+ 2*dt*(beta*V1(i)*X1(i)-k*Y1(i)*Z1(i)-delta*Y1(i)-F*Y1(i));
        Z2(i) = Z(i)+ 2*dt*(s*Y1(i)*Z1(i)-w*Z1(i)*Z1(i)-F*Z1(i));
        V2(i) = (4*V1(i)-V(i)+ 2*dt*delta*Y2(i))/(3+2*k0*dt*Z2(i)+2*gamma*dt);
        else
        X2(i) = X(i) - 2*dt*A*(X1(i+1)-X1(i-1))/2/dx+...
            2*dt*(lambda*X1(i)-beta*V1(i)*X1(i)-F*X1(i));
        Y2(i) = Y(i) - 2*dt*A*(Y1(i+1)-Y1(i-1))/2/dx+...
            2*dt*(beta*V1(i)*X1(i)-k*Y1(i)*Z1(i)-delta*Y1(i)-F*Y1(i));
        Z2(i) = Z(i) - 2*dt*A*(Z1(i+1)-Z1(i-1))/2/dx+...
            2*dt*(s*Y1(i)*Z1(i)-w*Z1(i)*Z1(i)-F*Z1(i));
        end
    end
    
    U2(1) = 0;
    for i = 1:N-1
        Fp = lambda*X2(i+1)-mu*(1-X2(i+1)-Y2(i+1)-Z2(i+1))+...
            s*Y2(i+1)*Z2(i+1)-w*Z2(i+1)^2;
        Fm = lambda*X2(i)-mu*(1-X2(i)-Y2(i)-Z2(i))+...
            s*Y2(i)*Z2(i)-w*Z2(i)^2;
        U2(i+1) = (xx(i)^2*U2(i)+R(iter+1)/2*dx*(xx(i+1)^2*Fp+xx(i)^2*Fm))/(xx(i+1)^2);
    end
    
    A1 = -(xx.*U2(N)/R(iter+1)+2*D./(R(iter+1)^2)./xx);
    A1(1) = 0;
    A2 = -D./R(iter+1)^2;
    S = 4*V1-V+2*delta*dt*Y2;
    L1 = 2*dt/dx/dx*A2-dt/dx*A1;
    D1 = 3-4*dt/dx^2*A2+2*k0*dt*Z2+2*gamma*dt;
    UU1 = 2*dt/dx^2*A2+dt/dx*A1;
    L1(N) = 4*dt/dx/dx*(A2)-dt/dx*A1(N);
    UU1(1) = 4*dt/dx^2*A2+dt/dx*A1(1);
    V2 = tridiagSolve(L1(2:N), D1(1:N), UU1(1:N-1), S);
    
    
    X = X1;Y = Y1;Z = Z1;U = U1;V = V1;
    X1 = X2;Y1 = Y2;Z1 = Z2;U1 = U2;V1 = V2;
 
    

    %V(:,iter+1) = tridiagSolve(L1, D1, U1, B1)
    if mod(iter,10) == 0
        X1 = 0.5*(X+X2);
        Y1 = 0.5*(Y+Y2);
        Z1 = 0.5*(Z+Z2);
        U1 = 0.5*(U+U2);
        V1 = 0.5*(V+V2);
    end
end
figure(2);plot(t,R);title('R')