function heat1d(J, N)

% initialization
    
    T = 0.2;
    t = 0.0;
    
    dx = 1/(J-1);
    dt = (T-t)/N;
    nu = dt/(dx*dx);
    
    x = linspace(0,1,J)';
    u = u_init(x);
    
    u1 = zeros(J,1);
    
% main loop for integration

    for n = 1:N
       for j = 2:J-1
          u1(j) = u(j) + nu*(u(j+1)-2*u(j)+u(j-1));
       end
       % for BCs
       u1(1) = 0.0;
       u1(J) = 0.0;
       % update time and solution
       t = t + dt;
       u = u1;
       % plot the solution
       figure(1)
       plot(x,u,'r-o');
       axis([0 1 0 1]);
       drawnow;
    end

% print parameters

    fprintf('\n%d %d %e %e %e\n', J, N, dt, dx, nu);
          
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = u_init(x)

    u = sin(pi*x);


    
    
    

