function CalorCondVar

L = 1; % Tamanio del dominio
N = 100; %Numero de sub-intervalos en la discretizacion
dx = L/N;
T = .5;
PlayMovie = 0;

%% Definir el eje x y conductividad

x = zeros(N+2,1); %Agregando las dos celdas fantasma
for j=0:N+1
    j_shift = j+1;
    x_j = dx*( j-1/2 );
    x(j_shift,1) = x_j;
end

sigma = .01;
nu_mhalf = zeros(N+2,1); nu_phalf = zeros(N+2,1);
for j=0:N+1
    j_shift = j+1;
    x_jphalf = dx*( j ); x_jmhalf = dx*( j-1 );
    nu_mhalf(j_shift,1) = x_jmhalf/L*0.1; % 0.1*exp(-(x_jmhalf-L/2)^2/sigma);
    nu_phalf(j_shift,1) = x_jphalf/L*0.1; %0.1*exp(-(x_jphalf-L/2)^2/sigma);
%     if L/2-L/8 < x_jmhalf && x_jmhalf < L/2+L/8
%         nu_mhalf(j_shift,1) = 0;
%     else
%         nu_mhalf(j_shift,1) = 0.1;
%     end
%     if L/2-L/8 < x_jphalf && x_jphalf < L/2+L/8
%         nu_phalf(j_shift,1) = 0;
%     else
%         nu_phalf(j_shift,1) = 0.1;
%     end
end
nu = (nu_mhalf+nu_phalf)/2;

nu_max = max(nu_mhalf);
nu_max = max(nu_max,max(nu_phalf));
nu_max = max(nu_max,max(nu));

% plot(x,nu)
% stop

%% Definir la condicion inicial

u0 = zeros(N+2,1);

for j=0:N+1
    j_shift = j+1;
    x_j = x(j_shift,1);
    if L/3 < x_j && x_j < 2*L/3
        u0(j_shift,1) = 1;
    else
        u0(j_shift,1) = 0;
    end
    %u0(j,1) = cos(sin(2*pi/L*x_j)^2);
end

% figure(1)
% plot(x,u0,'b -','LineWidth',2)
% xlabel('x','FontSize',20)
% ylabel('u_o','FontSize',20)
% title('Condicion inicial','FontSize',20)
% stop

%% Solucion numerica

u_n = u0; 
u_np1 = u_n;

R = 0.45;
dt = R*dx^2/nu_max;

t=0;
while t < T
    for j=1:N
        j_shift = j+1;
        nu_jmhalf = nu_mhalf(j_shift,1); 
        nu_jphalf = nu_phalf(j_shift,1); 
        nu_j = (nu_jmhalf+nu_jphalf)/2;
        u_np1(j_shift,1) = u_n(j_shift,1) + (dt*nu_jmhalf/dx^2)*u_n(j_shift-1,1)-2*(dt*nu_j/dx^2)*u_n(j_shift,1)+(dt*nu_jphalf/dx^2)*u_n(j_shift+1,1);
    end
    %Condicion de frontera en la izquierda: Dirchlet = 0
    j = 0; j_shift = j+1;
    u_np1(j_shift,1) = u_np1(j_shift+1,1);
    %Condicion de frontera en la derecha: Dirchlet = 0
    j = N+1; j_shift = j+1;
    u_np1(j_shift,1) = 0; %u_np1(j_shift-1,1);
    
    u_n = u_np1;
    t=t+dt
    
    if PlayMovie == 1
%         clf
%         figure(1)
%         plot(x(2:N+1),u_np1(2:N+1),'r -'); %,'LineWidth',2)
%         axis([0 L 0 1.5])
%         print('test','-depsc',figure(1))
        
        contornos = zeros(N+2,2);

        contornos(:,1) = u_np1;
        contornos(:,2) = u_np1;

        y = zeros(2,1);
        y(1,1) = -.1; y(2,1) = .1;

        figure(1)
        hold on
        contour(x(2:N+1,1),y,transpose(contornos(2:N+1,1:2)),100); colorbar;
        caxis([0 1])
        axis([0 L -1 1])
        print('test','-depsc',figure(1))
    end
    
end

%% Graficar condici?n inicial y soluci?n aproximada a tiempo T

figure(1)
hold on
plot(x(1:N+1),u_np1(1:N+1,1),'r --','LineWidth',2)
xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Soluciones numerica','FontSize',20)
lg = legend('Solucion numerica a tiempo T');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
axis([0 L 0 1.5])
%print('SolTpt01Rpt55','-depsc',figure(1))
stop

%% Contoour

contornos = zeros(N+2,2);

contornos(:,1) = u_np1;
contornos(:,2) = u_np1;

y = zeros(2,1);
y(1,1) = -.1; y(2,1) = .1;

figure(1)
hold on
contour(x(2:N+1,1),y,transpose(contornos(2:N+1,1:2)),100); colorbar;
caxis([-0.1 1.1])
axis([0 L -1 1])
stop

xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Soluciones numerica','FontSize',20)
lg = legend('Solucion numerica a tiempo T');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
axis([0 L 0 1.5])
%print('SolTpt01Rpt55','-depsc',figure(1))
