function DifusionDirichlet

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

%% Definir el eje x

x=zeros(N,1);
for j=1:N
    x_j = dx*( j-1/2 );
    x(j,1) = x_j;
end

%% Definir la condicion inicial

u0 = zeros(N,1);

for j=1:N
    x_j=x(j,1);
    if L/3 < x_j && x_j < 2*L/3
        u0(j,1) = 1;
    else
        u0(j,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 = zeros(N+2,1);

u_n(2:N+1) = u0; %Celdas dentro del dominio
u_n(1,1) = u_n(N+1,1); %Celdas fantasma con condiciones periodicas
u_n(N+2,1) = u_n(2,1);

u_np1 = u_n;

R = 0.49;
dt = R*dx^2/nu;

t=0;
while t < T
    for j=1:N
        j_shift = j+1;
        u_np1(j_shift,1) = u_n(j_shift,1) + dt*nu*(u_n(j_shift-1,1)-2*u_n(j_shift,1)+u_n(j_shift+1,1))/dx^2;
    end
    
    %Condiciones de frontera:
    %Izquierda: Dirichlet 1
    j = 0; j_shift = j+1;
    u_np1(j_shift,1) = 1;
    j = N+1; j_shift = j+1;
    %Izquierda: Dirichlet 2
    u_np1(j_shift,1) = 2;
    
    u_n = u_np1;
    t=t+dt
    
    if PlayMovie == 1
        clf
        figure(1)
        plot(x,u_np1(2:N+1),'r -.')
        print('test','-depsc',figure(1))
    end
end

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

figure(1)
hold on
plot(x,u_np1(2:N+1),'r --','LineWidth',2)
xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Soluciones exacta y numerica','FontSize',20)
lg = legend('Solucion a T = 10');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
print('SolDirichT10','-depsc',figure(1))
