clc;
clear all;

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Forward Euler caso DuDt  =  a*D2uDx2
% con u(0,t)=0
%     u(L,t)=0
%     u(x,0)=Tmax*4/L^2*(x*(L-x))
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
L  =  30;   % Longitud del intervalo espacial
T  =  300;  % Tiempo final de la simulación
Nx =  100;  % Number sub-intervals in space
Nt =  10000;  % Number sub-intervals in time

x    =  (linspace(0,L,Nx+1))';
time =  (linspace(0,T,Nt+1))';
h    =  L/Nx; % Espaciado para variable espacial x
Dt   =  T/Nt; % Espaciado para variable temporal t

a    =  1;  % Conductividad térmica del material

c    =  a*Dt/h^2;
fprintf('Courant number is %f\n',c)


u    =  zeros(Nx+1,Nt+1);
u0   =  0;
uL   =  0;
u(1,:)  =  u0;   % Condición de contorno en x=0
u(end,:)  =  uL;    % Condición de contorno en x=L

%u(:,1)    =  4/L^2*(x.*(L-x)); % Condición inicial
for ix=1:Nx+1
  u(ix,1) = 4/L^2*x(ix)*(L-x(ix));
end
%figure(1)
%plot(x,u(:,1),'LineWidth',4)
%xlabel('x')
%ylabel('u(x,0)')
%grid on

 

Matrix    =  eye(Nx-1) - c*MyToeplitz(Nx-1);
Vector    =  zeros(Nx-1,1);
Vector(1) =  u0;
Vector(end)  =  uL;
     

for itime=2:Nt+1
     u(2:end-1,itime)  =  Matrix*u(2:end-1,itime-1) + c*Vector;
end

%Snapshot  =  sort(randi(Nt+1,9,1));
%Snapshot(1)  =  1;
Snapshot     =  [1;1000;2000;3000;4000;5000;6000;7000;8000];
figure(1)
title('Forward Euler')
for iplot=1:9
subplot(3,3,iplot)
plot(x,u(:,Snapshot(iplot)),'LineWidth',4)
ylim([0 1]);
xlabel('x')
ylabel('u(x,t)')
grid on
end

tmesh   =  repmat(time,1,length(x));
xmesh   =  repmat(x',length(time),1);
umesh   =  reshape(u',Nt+1,Nx+1);

figure(2)
surf(xmesh,tmesh,umesh,'edgecolor','none')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
colorbar
  


%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Bacward Euler caso DuDt  =  a*D2uDx2
% con u(0,t)=0
%     u(L,t)=0
%     u(x,0)=Tmax*4/L^2*(x*(L-x))
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
L  =  30;
T  =  300;
Nx =  100;  % Number sub-intervals in space
Nt =  10000;  % Number sub-intervals in time

x    =  (linspace(0,L,Nx+1))';
time =  (linspace(0,T,Nt+1))';
h    =  L/Nx;
Dt   =  T/Nt;

a    =  1;

c    =  a*Dt/h^2;
fprintf('Courant number is %f\n',c)

u    =  zeros(Nx+1,Nt+1);
u0   =  0;
uL   =  0;
u(1,:)  =  u0;   % Condición de contorno en x=0
u(end,:)  =  uL;    % Condición de contorno en x=L
u(:,1)    =  (4)/L^2*(x.*(L-x)); % Condición inicial


Matrix    =  eye(Nx-1) + c*MyToeplitz(Nx-1);
Vector    =  zeros(Nx-1,1);
Vector(1) =  u0;
Vector(end)  =  uL;


for itime=2:Nt+1
    T                 =  u(2:end-1,itime-1) + c*Vector; 
    % Matrix*u=T
    u(2:end-1,itime)  =  Matrix\T; % Solver directo
     
end

%Snapshot  =  sort(randi(Nt+1,9,1));
%Snapshot(1)  =  1;
Snapshot     =  [1;1000;2000;3000;4000;5000;6000;7000;8000];

figure(3)
title('Backward Euler')
for iplot=1:9
subplot(3,3,iplot)
plot(x,u(:,Snapshot(iplot)),'LineWidth',4)
ylim([0 1]);
xlabel('x')
ylabel('u(x,t)')
grid on
end

tmesh   =  repmat(time,1,length(x));
xmesh   =  repmat(x',length(time),1);
umesh   =  reshape(u',Nt+1,Nx+1);

figure(4)
surf(xmesh,tmesh,umesh,'edgecolor','none')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
colorbar




