% This program is written by Chiu-Yen Kao
% for course 807
function u = heat_eqn()

nx = 11;
h = 1/(nx-1);
x = 0:h:1;
u = init_fcn(x);
ub0 = 0; % boundary condition
ub1 = 1; % boundary condition
figure(1);
plot(x,u,'r')
legend('initial condiion')
k = 0.0048;
nt = 0.1/k; % number of iteration
r = k/h^2
plot_iter = 0;
for iter = 1:nt
    t = iter*k;
    up = u(3:nx);
    um = u(1:nx-2);
    uc = u(2:nx-1);
    u(2:nx-1) = r*up+(1-2*r)*uc+r*um;
    if mod(iter,10)==0
        hold on;plot(x,u,'b');
    end
end

function val = init_fcn(x);
val = (2*x).*(x<=0.5)+(2*(1-x)).*(x>0.5);
