% yn=heateq(y0,t0,tn,N)
% 
% Solve the heat equation on the unit interval.
% For more information see Morton and Mayer
% Numerical Solution of Partial Differntial
% equations page 11.
%
%    y0 -- initial distribution of heat
%    t0 -- time of the initial condition
%    tn -- ending time
%    N  -- number of time steps to take

function yn=heateq(y0,t0,tn,N)
    J=length(y0)-1;
    dx=1/J;
    dt=(tn-t0)/N;
    mu=dt/(dx*dx);
    disp(sprintf('mu=%g',mu));
    yn=y0;
    for n=0:N-1
        ynp1(J+1)=0;
        ynp1(1)=0;
        for j=2:J
            ynp1(j)=yn(j)+...
            mu*(yn(j+1)-2*yn(j)+yn(j-1));
        end
        yn=ynp1;
    end