function BurgersViscNeumann

L = 3; % Tamanio del dominio
eps = .005*4; %Coeficiente de difusion
N = 1000; %Numero de sub-intervalos en la discretizacion
%dx = L/N;
T = 2; %1; %.5;

x_initial = -L;
x_final = L;
dx = (x_final - x_initial)/N;

nu = 0.45;
dt = nu*dx^2/eps;

%% Definir el eje x
x = zeros(N,1);
for j=1:N
    x_j = x_initial + 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 -1 <= x_j && x_j <= 0
        u0(j,1) = (-x_j-1);
    end
    if 0 <= x_j && x_j <= 1
        u0(j,1) = (-1+x_j);
    end
end

% plot(x,u0)
% stop

%% Aplicar la transformada de Cole-Hopf

% phi = exp( -1/(2*eps)*int_{0}^x u(s,t)ds )
phi = zeros(N,1);
integral = 0;
for j=1:N
    integral = integral + u0(j,1)*dx;
    phi(j,1) = exp(-1/(2*eps)*integral); 
end

% plot(x,phi,'r .')
% st


%% Solucion para la ecuacion del calor

phi_old = phi;
phi_new = phi;

t = 0; %cont=0;
while t < T
    for j=2:N-1
        phi_new(j,1) = phi_old(j,1) + eps*dt/dx^2*( phi_old(j-1,1) -2*phi_old(j,1) + phi_old(j+1,1) );
    end
    %Casos especiales
    %Condiciones de frontera Neumann
    %j=1 , celda 1
    phi_j0 =  phi_old(1,1);
    phi_new(1,1) = phi_old(1,1) + eps*dt/dx^2*( phi_j0 -2*phi_old(1,1) + phi_old(2,1) );
    %Celda j=N
    phi_jNp1 = phi_old(N,1);
    phi_new(N,1) = phi_old(N,1) + eps*dt/dx^2*( phi_old(N-1,1) -2*phi_old(N,1) + phi_jNp1 );
    
    phi_old = phi_new;
    
    t = t+dt;
    display(t)
    
%     if t >= T/10*cont
%         figure(1)
%         plot(x,phi_new)
%         print('Fig1','-depsc',figure(1))
%         cont = cont+1;
%     end
end
%stop

phi_new_x = zeros(N,1);
for j=2:N-1
    phi_new_x(j,1) = (phi_new(j+1,1)-phi_new(j-1,1))/(2*dx);
end
phi_new_x(1,1) = 0;
phi_new_x(N,1) = 0;

% plot(x,phi_new_x)
% stop

%% Volver a u con la inversa de Cole-Hopf

u_sol = -2*eps*phi_new_x./phi_new;

%% Solucion con diferencias finitas

u_old = u0;
u_new = u0;

max_evalues = max(abs(u0));
nu_Roe = 0.9;
dt = nu_Roe*dx/max_evalues;

t = 0;
while t < T
    
    for j=2:N-1
        lambda = (u_old(j,1)+u_old(j-1,1))/2;
        if lambda >= 0
            u_new(j,1) = u_new(j,1) - lambda*dt/dx*( u_old(j,1) - u_old(j-1,1) );
        end
        lambda = (u_old(j,1)+u_old(j+1,1))/2;
        if lambda <= 0
            u_new(j,1) = u_new(j,1) - lambda*dt/dx*( u_old(j+1,1) - u_old(j,1) );
        end
    end
    u_new(1,1) = 0;
    u_new(N,1) = 0;
    
    u_old = u_new;
    
    t = t+dt;
    display(t)
end

% plot(x(2:N-1,1),u_sol(2:N-1,1),'b -','LineWidth',2)
% 
% stop

%% Graficar condicion inicial y solucion aproximada a tiempo t

figure(1)
hold on
plot(x(2:N-1,1),u0(2:N-1,1),'b -','LineWidth',2)
plot(x(2:N-1,1),u_sol(2:N-1,1),'r -','LineWidth',2)
plot(x,u_new,'r --','LineWidth',2)
xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Condicion inicial y solucion','FontSize',20)
lg = legend('Condicion inicial','Solucion Viscosa a tiempo T','Soluci?n referencia');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
print('FigT2Viscpt02','-depsc',figure(1))

