clear all
close all
%------------------------------
% Solution of the 1d heat equation
% u_t=D*u_xx
% D-diffusion coefficient;
% Initial condition: u(0,x)=f(x)-itinial heat distribution;
% x in (-L,L) with
% no-flux boundary conditions: u_x(t,-L)=0, u_x(t,L)=0.
%-----------------------------
% The scheme: the BTCS implicit method 
% u^n+1=A^-1*b
% a=D*dt/dx^2;
%-----------------------------
%-----------------------------
L=5; % length of the string
D=1;
M=100; % amount of points on the (0,L)
x=linspace(-L, L, M+1);
dx=x(2)-x(1);% space grid size
dt=0.05; % time grid size
a=(D*dt)/dx^2;
T=200; % amount of time steps
tt=20;
ys=1.1;
rad=0.01;
nn=1;
a1=-a;
a2=1+2*a;
a3=-a;
b=exp(-L^2/4)/sqrt(4*pi);
%-------------------------
v=zeros(M+1,1);
% v(1)=-2*a*dx*b;
% v(M+1)=-2*a*dx*b;
%-----------------------------
u=zeros(T+1,M+1);


f=inline('exp(-y.^2)','y'); % an initial wave u(x,)

u(1,:)=f(x);

plot(x,u(1,:),'r','LineWidth',2)
hold on
%plot(x,-1.5*x+2,'r','LineWidth',2)
axis([-L-rad L+rad 0 ys])
       title(['t=' num2str(1)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       pause(nn)
       hold off
%-----------------------------
%-----------------------------
% Calculate the matrix A, size(A)=M-1xM-1
ah=a2*ones(M+1,1);
ao=a3*ones(M,1);
au=a1*ones(M,1);
A=diag(ah)+diag(au,-1)+diag(ao,1);
A(1,2)=-2*a;
A(M+1,M)=-2*a;
%-----------------------------
%Calculate all other time rows
for t=1:T
    pru=u(t,:)'+v;
    u(t+1,:)=inv(A)*pru;
    if mod(t,tt)==0
       plot(x,u(t+1,:),'r','LineWidth',2);
       hold on
       %plot(x,exp(-x.^2/(4*pi*(t+0)*dt))/sqrt(4*pi*(t+0)*dt),'r','LineWidth',2);
       axis([-L-rad L+rad 0 ys])
       title(['t=' num2str(t)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       hold on 
       pause(nn)
       hold off
    end;
end;
tz=0:T;
figure(2)
[C,h]=contourf(x,tz,u,100);
axis square
set(h,'LineStyle','none')
xlabel('x')
ylabel('t')
colorbar
grid off

figure(3)
h=surf(x,tz,u)
shading interp
light('Position',[-2,2,20])
lighting phong
material([0.4,0.6,0.5,30])
set(h,'EdgeColor','none','FaceColor',[0.2,0.2,0.9],'BackFaceLighting','lit')
alpha(0.8)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
axis([-L-rad L+rad 0 T 0 1])
