heat.m

function heat(t_i,t_f,a,l,dx,dt,c)
    n      = (t_f-t_i)/dt;
    m      = (l-a)/dx;
    lambda = c^2*dt/dx^2;
    u      = zeros(m+1,n+1);

    x = a:dx:l;
    t = t_i:dt:t_f;

    for j=1:m+1,
        u(j,1) = sin(x(j));     % u(x,0) = g(x)
    end

    for k=1:n
        u(1,k+1)=0.;            % u(0,t) = 0
        for j=2:m
            u(j,k+1)=u(j,k)+lambda*(u(j+1,k)-2*u(j,k)+u(j-1,k))+dt*x(j);
        end
        u(m+1,k+1)=0.;          % u(l,t) = 0
    end

    % exact solution comparison
    for j=1:m
        for k=1:n
            [j k x(j) t(k) u(j,k) exp(-3*t(k))*sin(x(j))+x(j)*t(k)]
        end
    end