%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% HeatEquation.m
%
% This simple script outputs plots of numerical solutions to the heat equation
% using a simple semi-discretization with centered finite
% differences. The initial profile must be changed manually by the user.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Subroutine to collect initial input
%n=the number of spatial nodes. Input from user.
%k=diffusion constant. Input from user.

n=100;
k=.01;
kappa=1;
r=.01;
tmax=100;
timepoints=100;

%% Construction of spatial domain and spacing
%x=the spatial domain.
x=linspace(-4,4,n);

%dx=spacing of nodes
dx=x(2)-x(1);


%u=The initial profile of the wave. We set the profile to be a Gaussian function

u0=2*sin(x).^2;

%tspace= is the times we want to capture the wave as it moves. Note that we
%are using a very course spacing but the Matlab ode solver will
%automatically generate a finer temporal mesh in its routines.
tspace=linspace(0,tmax,timepoints);

%% Numerically Solving PDE
%These options tell me how long the solver runs for.
options=odeset('Stats','on');

%Construction of second order spatial derivative.
D2=diag(-2*ones(1,n),0)+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1);
D2=D2/dx^2;

%Enforcement of boundary conditions
D2(:,1)=0;
D2(:,n)=0;



%Numerically solve heat equation. 
[t,u]=ode45(@ReactionDiffusion_diff,tspace,u0,options,k,kappa,r,D2);

%y is a matrix. The i row corresponds to the ith time we are 

%Plot the solution at the selected times.
for i=1:timepoints
    u(i,1)=0; %Rigid Enforcement of B.C.
    u(i,n)=0;
    plot(x,u(i,1:n),'k-');
    axis([-4,4,0,2]);
    temp=strcat('Diffusion Equation, Gaussian, t=',num2str(t(i)));
    title(temp);
    xlabel('x');
    ylabel('u');
    drawnow;
end


