function EcuacionDifusion

L = 1; % Tama?o del dominio
N = 1000; %N?mero de sub-intervalos en la discretizaci?n
dx = L/N;
T = 100;
%Conductividad: e^(-alpha*t)
alpha = 100;

%% Definir el eje x

x=zeros(N,1);
for j=1:N
    x_j = dx*( j-1 );
    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

% plot(x,u0)
% stop


%% Descomposici?n 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

% plot(an(1:500,1),an(1:500,2))
% stop

%% 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(-(2*pi/L*n)^2/alpha*(1-exp(-alpha*T)));
    end
end

% u0comp = real(u0comp);
% 
% plot(x,u0comp)
% stop

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

figure(1)
hold on
plot(x,u0,'b -','LineWidth',2)
plot(x,u_sol,'r -','LineWidth',2)
xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Condicion inicial y solucion','FontSize',20)
lg = legend('Condicion inicial','Solucion a tiempo T = 100');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
print('SolAlpha100T100','-depsc',figure(1))
