
%                          Matlab program heat.m
%
%   This program integrates the heat equation  u_t -ku_xx = 0
%   on the interval [0,10] using the Crank-Nicholson method.
%   The implicit set of equations are solved at each time step
%   using an LU factorization. The time step delt = .05, and 
%   delx = .05 as well.



k = 1; % the value of the diffusion coefficient k 

f = @(x) sin(x); % initial condition

left = @(x) 0*x; % left boundary condition
right = @(x) 0*x % right boundary condition


delt = .05
disp(' Enter the final time for the computation ')
T = input (' T = ')

N=T/delt;

Nframes=200;
heatmovie=moviein(Nframes-1);

nframe=floor(N/Nframes);
N=nframe*Nframes;

delx = .05
r = .5*k*delt/(delx^2);

J= 10/delx;

x = 0:delx:10;

    snap0 = f(x);

v = snap0(2:1:J);

D = sparse(1:J-1, 1:J-1,(1+2*r), J-1, J-1, J-1);
E = sparse(2:J-1, 1:J-2, -r, J-1, J-1, J-1);
A = D + E +E';

[L,U] = lu(A);


% Note: the original index j runs from j = 1 ( x = 0) to j = J ( x = 10 -delx).
% The index in v (and hence in b and z) runs from j = 1 (x= delx) to
% j = J-1 (x = 10-delx).

for j=1:Nframes;
for n = 1:nframe
    b = 2*v ;
    b(1) = b(1) + r*( left(n*delt) +left((n-1)*delt) );
    b(J-1) = b(J-1) + r*( right(n*delt) +right((n-1)*delt) ) ;
    y = L\b';
    z = (U\y)';
    v = z- v;
end 

snap1 = [left(n*delt),v,right(n*delt)];
disp('Computed at time t='),n*j*delt


plot(x,snap0,x,snap1), hold off
movie(:,j)=getframe;
end