%% Euler forward Dirichlet
L = 1;  % domain size (in the lecture notes this is denote h)
T = 2;  % time limit (max time)
f = @(x,t) 0*x.*t;    % rhs of the more general equation `u_t - d u_xx = f`
c1 = @(t) 1+0*t;      % _right_ boundary condition
c2 = @(t) 0*t;        % _left_ boundary condition
u0 = @(x) 0*x;        % initial values
D = 0.5;              % diffusion parameter `d` in `u_t - d u_xx = f`
%uex = @(x,t) cos(x).*exp(t);

N = 10;               % nr. of _space_ discretization points
K = 200;              % nr. of _time_ discretization points
[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);

% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
Delta_T = T / K;
Delta_x = L / N;
d = D * Delta_T / Delta_x^2;
fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
if d > 0.5
    fprintf("-> NOT Stable\n")
else
    fprintf("-> Stable\n")
end

figure(1)
for ii = 1:K+1 % iterates time
    hold on
    plot(x, u(:, ii)');
    xlim([0 L])
    pause(0.05);
    hold off
end

% 3D plot of space solution over time
space = linspace(0,L,101);
time = linspace(0,T,201);
[xx,yy] = meshgrid(time,space);
%exsol = uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)

%% Eulero forward Mixed BC

L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1.1;
%uex=@(x,t) cos(x).*exp(t);

N=25;
K=200;
[x,t,u]=Mixed_EA(L,N,T,K,c1,c2,f,u0,D);

figure(1)
for ii=1:K+1 
    plot(x,u(:,ii)');
    xlim([0 L])
    ylim([0 1.5])
    pause(0.02);
end

space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)

%% Eulero Backward Mixed BC

L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1;
%uex=@(x,t) cos(x).*exp(t);

N=25;
K=200;
[x,t,u]=Mixed_EI(L,N,T,K,c1,c2,f,u0,D);

figure(1)
for ii=1:K+1 
    plot(x,u(:,ii)');
    xlim([0 L])
    ylim([0 1.5])
    pause(0.02);
end

space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)