clear all
close all
clc

set(0,'DefaultLineLineWidth',2);
set(0,'DefaultAxesFontSize',16);
set(0,'DefaultLineMarkerSize',10);

%%%%%%%%%%%%% Physical parameters
L = 1;
Phi_w =  0;   % Boundary condition
Phi_0 = 40;   % Initial condition
Gamma = 0.1;
rho = 1;
t_End = 5;

%%%%%%%%%%%%% Numerical parameters
nx = 51;

% Specify number of time steps OR size of time step
%nt = 900;    dt = t_End/nt,
 dt = 0.02,   nt = t_End/dt,

theta = 0.5;   % for ex. 0 (explicit), 0.5 (Crank-Nicolson) or 1 (implicit)

%%%%%%%%%%%%%% Grid generation
dx = L/(nx-1);
Dx = dx;

x0 = 0:dx:L;
time = 0:dt:t_End;

check_dt_max_diffusion = rho*dx^2/Gamma,

%%%%%%%%%%%%%% Initial condition
Phi_init = Phi_0*ones(nx,1); 

%%%%%%%%%%%%%% Initialize the matrix
A = zeros(nx,nx);
b = zeros(nx,1);
for i=2:nx-1
    A(i,i-1) =              -theta*Gamma/dx; % W
    A(i,i+1) =              -theta*Gamma/dx; % E
    A(i,i)   = rho*dx/dt + 2*theta*Gamma/dx; % P
end

% Boundary conditions
A(1,1)   = 1;
b(1)     = Phi_w;
A(nx,nx) = 1;
b(nx)    = Phi_w;

%%%%%%%%%%%%%% Time marching
nt = length(time);
Phi = zeros(nx,nt);
Phi(:,1) = Phi_init;

col = hot(12);
j = 1;
figure('color','w'), plot(x0, Phi(:,1), 'color',col(j,:)),
hold on, grid on, box on,
xlabel('x'), ylabel('phi'), ylim([0 45]), drawnow, %pause
for k=2:nt    
    for i=2:nx-1        
        b(i) = Phi(i+1,k-1) * (1-theta)*Gamma/dx...
             + Phi(i-1,k-1) * (1-theta)*Gamma/dx...
             + Phi(i,  k-1) * (rho*dx/dt -2*(1-theta)*Gamma/dx);
    end
        
    % Solve
    Phi(:,k) = A\b;  
    
    if (mod(k, round(nt/5))==0)
        %hold off, 
        j = j+1;
        plot(x0, Phi(:,k), 'color',col(j,:)), hold on, grid on, box on,
        xlabel('x'), ylabel('\phi'), ylim([0 40]), drawnow
    end
end
exportgraphics(gca,'pure_diffusion-snapshots.pdf')


%%%%%%%%%%%%%% Plot solution T(x,t)
%figure('color','w','position',[874   378   560   420]), 
figure, hold on
xlabel('t'), ylabel('x'),
surf(time,x0, Phi, 'edgecolor','none')
c = colorbar;   c.Label.String = '\phi';
drawnow
exportgraphics(gca,'pure_diffusion-space_time_evolution.pdf')
 

%%%%%%%%%%%%%% Plot solution T(t) at x=L/2
figure('color','w'), hold on, grid on, box on
plot(time, Phi((nx+1)/2,:), 'b.')

%--- Theoretical solution (for constant initial condition)
time_theo = linspace(0,t_End,1001);
n_sum_theo = 500;
Sn = zeros(1, length(time_theo));
xp = L/2;
for j=1:n_sum_theo
    Sn = Sn + (1-(-1)^j)/j * sin(j*pi/L*xp) * exp(-Gamma/rho*(j*pi/L)^2*time_theo);    
end
Phi_theo = Phi_w + 2*(Phi_0-Phi_w)/pi*Sn;
plot(time_theo, Phi_theo, 'r--')

xlabel('t'), ylabel(['\phi(x=' num2str(x0((nx+1)/2)) ')'])
legend('Numerical','Theoretical')
exportgraphics(gca,'pure_diffusion-time_evolution.pdf')

title(['L=1, Phi_0=' num2str(Phi_0) ', Phi_w=' num2str(Phi_w)...
    ', \Gamma=' num2str(Gamma) ', \rho=' num2str(rho) ', n=' num2str(nx)...
    ', \Deltat=' num2str(dt) ', \theta=' num2str(theta)])