clear all
close all
%------------------------------
% Solution of the 1d heat equation
% u_t=D*u_xx=0
% c-diffusion coefficient;
% Initial condition: u(0,x)=f(x)-itinial heat distribution;
% x in (0,L) with
% boundary conditions: u(t,0)=Tl, u(t,L)=Tr.
%-----------------------------
% The scheme: explicit FTCS method 
% u_ij+1=au_i+1^j+(1-2a)u_i^j+au_i-1^j
% a=D*dt/dx^2;
%-----------------------------
%-----------------------------
L=1; % length of the string
D=1;
M=10; % amount of points on the (0,L)
x=linspace(0, L, M+1);
dx=x(2)-x(1);% space grid size
a=0.5;
dt=(a*dx^2)/D; % time grid size
T=40; % amount of time steps
tt=2;
ys=1.1;
Tl=0;
Tr=0;
rad=0.01
nn=1
%-----------------------------
u=zeros(T+1,M+1);

f=inline('4*y-4*y.^2','y'); % an initial wave u(x,)
%f=inline('-1.5*y+2+sin(pi*y)','y'); % an initial wave u(x,)

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

u(1,1)=Tl; u(1,M+1)=Tr;
plot(x,u(1,:),'k','LineWidth',2)
hold on
axis([0-rad L+rad 0 ys])
       title(['t=' num2str(1)],'Color','b')
       xlabel('x')
       ylabel('u(x,t)')
       pause(nn)
       hold off
%-----------------------------
% %Calculate all other time rows
for t=1:T
    u(t+1,1)=Tl; u(t+1,M+1)=Tr;
    for i=2:M
        u(t+1,i)=a*u(t,i-1)+(1-2*a)*u(t,i)+a*u(t,i+1);
    end;
    if mod(t,tt)==0
       plot(x,u(t+1,:),'r','LineWidth',2);
       hold on
       axis([0-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;

% y=0:0.01:1;
%       figure(2)
%       plot(y,(pi*y).^2 ,'k','LineWidth',2 )
%       hold on
%       plot(y,-log(1-4*0.1*sin((pi*y)/2).^2)/0.1 ,'r','LineWidth',2 )
%       hold on
%       plot(y,-log(1-4*0.2*sin((pi*y)/2).^2)/0.2 ,'m.','LineWidth',2 )
%       hold on
%       plot(y,-log(1-4*0.3*sin((pi*y)/2).^2)/0.3 ,'b','LineWidth',2 )
%       hold on
%       plot(y,-log(1-4*0.4*sin((pi*y)/2).^2)/0.4 ,'g-.','LineWidth',2)
%       xlabel('k\Delta x/\pi')
%       ylabel('i\omega\Delta t/\alpha')
%      legend('exact','\alpha=0.1','\alpha=0.2','\alpha=0.3','\alpha=0.4')

% tz=0:T;
% figure(2)
% [C,h]=contourf(x,tz,u,100);
% set(h,'LineStyle','none')
% xlabel('x')
% ylabel('t')
% grid off
% colorbar

% tz=1:tt:T+1;
% figure(3)
% h=surf(x,tz,u(tz,:))
% 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.9)
% xlabel('x')
% ylabel('t')
