  function [M1, M2, F, b, x, e_conn] =  assemb_co(param)
  % The FEM equation for the temp. dist at time t_{n+1} satisfies
  %   (M_1 + M_2) z^{n+1} - M_1 z^n + F + b = 0 
    
%% Initialization & geometry
%----lines  omitted
%% Set up codistributed structure

% column pointers and such for codistributed arrays
  Vc  = codcolon(1, n_equations);
  lP  = localPart(Vc); lP_1 = lP(1); lP_end = lP(end);
  dPM = distributionPartition(codistributor(Vc));
  col_shft = [0 cumsum(dPM(1:end-1))];
  
% local sparse arrays
  M1_lab = sparse(n_equations, dPM(labindex));  M2_lab = M1_lab;
  b_lab  = sparse(dPM(labindex), 1); F_lab = b_lab;

%% Build the finite element matrices - Begin loop over elements
  for n_el=1:n_elements
    nodes_local     = e_conn(n_el,:);% which nodes are in this element
    % subset of nodes/columns on this lab
    lab_nodes_local = my_extract( nodes_local, lP_1, lP_end);
    if ~isempty(lab_nodes_local) % continue the calculation for this elmnt
%--- calculate local arrays - lines omitted