function EcuacionDifusion

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

%% 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

%% Descomposicion en modos de Fourier

M = N/2; %N?mero de modos de Fourier. La suma va de -M a M

%Calcular los a_n's
an = zeros(2*M+1,2);

%La primer entrada tendr? a_o
% De la 2 a la M+1 tendr? de a_1 a a_M
% De la M+2 a la M+2+M-1=2M+1 lleva de a_{-M} a a_{-1} 
i=sqrt(-1);
for n = -M:M
    
    % Calcular $ (1/L)*\int_0^L u_o(x)*exp(-2*\pi*i/L*n*x) dx $
    a_n = 0;
    for j=1:N %Suma de Riemann en cada celda
        x_j = x(j,1);
        a_n = a_n + (1/L)*u0(j,1)*exp(-2*pi*i/L*n*x_j)*dx;
    end
    
    if n >= 0 
        emf = n+1;
    else
        emf = 2*M+1 + n+1;
    end
    an(emf,1) = n;
    an(emf,2) = a_n;
end

%% Composici?n de los modos de Fourier y la soluci?n a la ecuacion de difusi?n

u0comp = zeros(N,1);
u_sol = zeros(N,1);

for emf = 1:2*M+1
    
    n = an(emf,1);
    a_n = an(emf,2);
    
    for j=1:N
        x_j = x(j,1);
        u0comp(j,1) = u0comp(j,1) + a_n*exp(2*pi*i/L*n*x_j);
        u_sol(j,1) = u_sol(j,1) + a_n*exp(2*pi*i/L*n*x_j)*exp(-nu*(2*pi/L*n)^2*T);
    end
end

%u0comp = real(u0comp);

%% 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
    j = 0; j_shift = j+1;
    u_np1(j_shift,1) = u_n(j_shift,1) + dt*nu*(u_n(N+1,1)-2*u_n(j_shift,1)+u_n(j_shift+1,1))/dx^2;
    j = N+1; 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(2,1))/dx^2;
    
    u_n = u_np1;
    t=t+dt
end

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

figure(1)
hold on
plot(x,u_sol,'b -','LineWidth',2)
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 exacta','Solucion numerica a tiempo T');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
%print('SolTpt01Rpt55','-depsc',figure(1))
