% Dirchlet Boundary Conditions

%% Parameters:

a=0; %Spatial Domain
b=1;

n=100; %Number of spatial points
m=10000; %Number of temporal points

gamma=1; %Diffusion constant

%% Constuction of spatial-temporal mesh

x=linspace(a,b,n);
dx=x(2)-x(1);

dt=1/10*(dx^2/gamma); %Satisfying CFL condition
t=linspace(0,(m-1)*dt,m);

%% Initial profile
u0=sin(2*pi*(x-a)/(b-a)).^3.*sin(sqrt(2)*x);

%% Construction of second derivative matrix

upperdiag=ones(1,n-1);
lowerdiag=ones(1,n-1);
centerdiag=-2*ones(1,n);

D=diag(upperdiag,1)+diag(centerdiag,0)+diag(lowerdiag,-1);


D=1/(dx^2)*D;

D=sparse(D);

%% Method of lines

u=zeros(n,m); %Initialization of distribution
u(:,1)=u0;

vidObj = VideoWriter('Dirichlet.avi');
open(vidObj);

%Plotting Initial guess
figure('units', 'normalized', 'outerposition', [0,0,1,1]);
hold on;
plot(x,u0,'LineWidth',2);
ylabel('$u(x,t)$','Interpreter','latex','Fontsize',18);
xlabel('$x$','Interpreter','latex','Fontsize',18);
title(strcat('$t=$',num2str(t(1))),'Interpreter','latex','Fontsize',18 );
axis([a,b,-1,1]);
drawnow;

currFrame = getframe(gcf);
writeVideo(vidObj,currFrame);

% Looping through Changes in time
for i=1:m-1
    u(:,i+1)=u(:,i)+(gamma*D*u(:,i))*dt;
    
    u(1,i+1)=0; %boundary conditions
    u(n,i+1)=0; %boundary conditions
    
    if mod(i,100)==1, %Plot every 100th iterate
        clf;
        plot(x, u(:,i+1),'LineWidth',2);
        ylabel('$u(x,t)$','Interpreter','latex','Fontsize',18);
        xlabel('$x$','Interpreter','latex','Fontsize',18);
        title(strcat('$t=$',num2str(t(i))),'Interpreter','latex','Fontsize',18 );
        axis([a,b,-1,1]);
        drawnow;
        
        currFrame = getframe(gcf);
        writeVideo(vidObj,currFrame);
    end
    
end

close(vidObj);
