clc;
close all;
clear variables;

%geometry -  xi direction
x0=0;%[m]
x1=1;%[m]
nx=50;
x=linspace(x0,x1,nx);
dx=x(2)-x(1);

%material properties
lambda_x=2.2;%[W/mK]
cp=880;%[J/kgK]
rho=2200;%[kg/m3]
a=lambda_x/(cp*rho);%[m2/s]

%time
t_length=60*60*24;%[s]

%initial condition
T_0=20;%[°C]

%boundary conditions
T_E=0;%[°C]
T_I=20;%[°C]

%numerical settings
scheme='explicit';

%---------------------------------------------------------------------
switch scheme
    case 'explicit'
        %stable timestep
        dt_max=0.5*dx^2/a;
        nt=ceil(t_length/dt_max);
        t=linspace(0,t_length,nt);
        dt=t(2)-t(1);
        
        T=zeros(nx,nt);
        T(2:end-1,1)=T_0;
        T(1,1)=T_E;
        T(end,1)=T_I;
        
        %K matrix
        K=zeros(nx,nx);
        K(1,1)=1;
        for i=2:nx-1
            K(i,i-1)=a*dt/dx^2;
            K(i,i)=1-2*a*dt/dx^2;
            K(i,i+1)=a*dt/dx^2;
        end
        K(end,end)=1;
        %spy(K);
        
        %calculation
        for j=2:nt
           T(:,j)=K*T(:,j-1); 
        end
        
    case 'implicit'   
end

figure;
hold on;
for j=1:10:nt
    plot(x,T(:,j));
end


figure;
plot(t,T(5,:));

figure;
[X,Y]=meshgrid(x,t);
surf(X,Y,T','FaceColor','interp');