function temple3044_heateqn
%TEMPLE3044_HEATEQN
%    Solves the heat equation u_t=u_xx by finite differences.
%    The example uses forward Euler in time, central
%    differences in space, and Dirichlet boundary conditions.
%    The update rule is implemented using sparse matrices,
%    so the time stepping can easily be modified to implicit
%    methods.

% 03/2012 by Benjamin Seibold
%           http://www.math.temple.edu/~seibold/
% Feel free to modify for teaching and learning.

n = 25; % number of spacial grid cells
dt = 8e-4; % time step

n = 30;
%dt = .5/n^2
dt = 5.6e-4;

tf = 1; % final time
h = 1/n; % spacial resolution
nt = ceil(tf/dt); % number of time steps
dt = tf/nt; % actual time step
x = (0:h:1)'; % spacial grid
u0 = ic(x); % initial conditions
u = u0(2:end-1); % interior points only
A = spdiags(ones(n-1,1)*[1 -2 1]/h^2,-1:1,n-1,n-1); % Laplace matrix
I = speye(n-1); % identity matrix
b = zeros(n-1,1); % initialization for boundary condition vector

% Time loop
for j = 1:nt
    t = (j-1)*dt; % current time
    ub = bc(t); % boundary values at current time
    b([1 end]) = ub/h^2; % assign b.c. into r.h.s. vector
    u = (I+dt*A)*u+dt*b; % foward Euler time step
    % Plot solution
    clf, plot(x,u0,'b:',x,[ub(1);u;ub(end)],'r.-')
    title(sprintf('Solution at t=%0.2f',t+dt))
    xlabel('x'), ylabel('u'), drawnow
end

%==========================================================
function y = ic(x)
% Initial condition function
y = (1-(2*x-1).^2).^2;

function ub = bc(t)
% Boundary condition function
ub = [0,sin(4*pi*t).^2];
