function heattransfer_model(image_data, diffusivity_x, diffusivity_y, time_steps, output_file_path, divisor)
    % Turn image into grayscale array
    [X,map] = rgb2ind(image_data, 256);
    T = double(ind2gray(X,map));
    % Take a mean of the image data.
    %%T = mean(image_data, 3);
    flat_size = size(T);
    width = flat_size(1);
    height = flat_size(2);
    
    % Create two tempature array's
    T1 = T;
    T2 = zeros(width, height);
   
    %% initialize video output
    fig1 = figure(1);
    p = get(gcf, 'Position');
    set(gcf, 'Position', [p(1) p(2) height width] );
    set(gca, 'Position', [0 0 1 1] );
    winsize = get(fig1, 'Position');
    winsize(1:2) = [0 0];
    ht = moviein(time_steps/divisor, fig1, winsize);
    set(fig1, 'NextPlot', 'replacechildren');
    
    % Loop through each step in time
    for v=1:time_steps
        
        for i=2:(width-1)
            for j=2:(height-1)
                % Calculate new T2 from T1 for each mode
                T2(i,j) = T1(i,j) + diffusivity_x*(T1(i+1,j) + T1(i-1,j) - 2*T1(i,j)) + diffusivity_y*(T1(i,j+1) + T1(i,j-1) - 2*T1(i,j));
            end
        end
        T1 = T2;
        if mod(v,divisor) == 0
            image(T1), colormap(jet(256));        
            ht(:,v) = getframe(fig1, winsize);
        end
    end
    
    % Close video file?
    movie2avi(ht, output_file_path); 
end