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

NN = 200
N = 2*NN+1; % number of point in x
dx = 1/(N-1); 
xx = (0:(N-1))'*dx;
simp_coeff = ones(1,N)*dx/3;
simp_coeff(3:2:N-2) = 2*dx/3;
simp_coeff(2:2:N-1) = 4*dx/3;
simp_coeff = simp_coeff';

dt = 0.5*dx;%0.5*(dx)^2;
max_iter = N*2*10;
t = (0:(max_iter-1))*dt;
D =1;
R = zeros(1,max_iter);
V = zeros(N,1);
V1 = zeros(N,1);
V2 = zeros(N,1);
mu = 1;
gamma = 1;

% initialization


Vbar = 0.5;
Vtilde = 0.3;
Rs = fzero(@(x) tanh(x)-x/(1+((Vtilde/3/Vbar)*x^2)), 10);

R(1) = Rs; % mm at t=0
Vs = Vbar*Rs/sinh(Rs)*sinh(Rs*xx)./(Rs*xx);
Vs(1) = Vbar*Rs/sinh(Rs);

plot(Vs)

pause
%V = Vs;

V = 0.5*ones(N,1);


 
% first iteration
iter = 1;
a = 2.1812;



rat = R(iter)*sum((mu*(V-Vtilde).*xx.^2).*simp_coeff);
R(iter+1) = R(iter)+dt*rat;


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


%plot(abs(V1 - Vs))

for iter = 2:max_iter-1
    if mod(iter,20)==0
    figure(1);plot(xx,V1)
    legend('V')
    title(['T = ' num2str((iter-1)*dt) '    R = ' num2str(R(iter))]);
     drawnow
    end
    rat = R(iter)*sum(mu*(V-Vtilde).*xx.^2.*simp_coeff)
    R(iter+1) = R(iter)+dt*rat;

    A1 = -(xx.*rat/R(iter+1)+2*D./(R(iter+1)^2)./xx);
    A1(1) = 0;
    A2 = -D./R(iter+1)^2;
    S = 4*V1-V;
    L1 = 2*dt/dx/dx*A2-dt/dx*A1;
    D1 = (3-4*dt/dx^2*A2+2*gamma*dt)*ones(1,N);
    UU1 = 2*dt/dx^2*A2+dt/dx*A1;
    UU1(1) = 4*dt/dx^2*A2+dt/dx*A1(1);
    S(N-1) = S(N-1)-A1(N-1)*2*dt/2/dx*Vbar-A2*2*dt/dx/dx*Vbar;
    V2(1:N-1) = tridiagSolve(L1(2:N-1), D1(1:N-1), UU1(1:N-2), S(1:N-1));
    V2(N) = Vbar;
    
    V = V1;
    V1 = V2;
 
    
    %if mod(iter,10) == 0
    %    V1 = 0.5*(V+V2);
    %end
end
figure(2);plot(t,R);title('R'); axis([0 max(t) min(R)-1 max(R)+1])