JVAL=4; % JVAL is the J value of the SMBE subgraph
Ku=zeros(JVAL, JVAL+1); Km=zeros (JVAL, JVAL+1);
Ki=zeros(JVAL, JVAL+1); Kr=zeros (JVAL, JVAL+1);
Kc=zeros (JVAL, JVAL+1);
for n=1:JVAL, %n is the row index for the K matrices
    % Part I: compute the residue basis set
    a=1:n;
    a=a.*(2*a-1)% a is the vector of 2j choose 2 kappa coefficients
    RES=[]; % initialize RES as an empty matrix
    for j=n:-1:1, % start at the back of the flow graph and come in
        A=poly(-a(j:n)) % convert root form to power form
        [r p]=residue(1, A); % get the residues, 1=numerator, A=denom
        res=[zeros(1, j) fliplr(r')];
        RES=[RES; res]
    end %Note that the products of kappa’s cancel out
    A=poly(-[0 a]); % now do the last row in RES, i.e., j=0
    [r p]=residue (1, A);
    res=[fliplr(r')];
    RES=[RES; res];
    % Part II: use the residue basis set and SMBE graph to find
    % u_a(t) and m(t) given that the site contains exactly n actives.
    UT=n*RES(1,:); MT=0*res; % initialize the running totals
    IT=0*res; %incomplete
    RT=0*res; %reciprocal
    CT=0*res; %complex
    M=[zeros(1, n) 1]; % start with 1 in the nth node
    for m=1:(n+1),
        nM=zeros (m+1, n+1); % stack at the next m value
        for J=(n+1):-1:2,
            I=1;
            while (I<=m)&(I<=J),
                i=I-1; % i and j are the math indices
                j=J-1; % I and J are the matrix indices
                nm1=(j-i); % (eurejoin: eu within hopefuls)
                nm2=2*(j-i)*(j-i-1); %(mis within hopefuls)
                nm3=i*(2*i-1); %(mis within destined)
                nm4=4*(j-i)*i; %(mis between hopefuls and destineds)
                M(I, J-1)=M(I, J-1)+nm1*M(I, J);
                UT=UT+(j-1)*M(I, J-1)*RES(n+2-j,:);
                MT=MT+(m-1)*M(I, J-1)*RES(n+2-j,:);
                if ((m-1)==1) IT=IT+(m-1)*M(I, J-1)*RES(n+2-j,:);
                elseif ((m-1)==2)&((j-1)<2)&(i==0)
                    RT=RT+(m-1)*M(I, J-1)*RES(n+2-j,:);
                else CT=CT+(m-1)*M(I, J-1)*RES(n+2-j,:);
                end
                nM(I, J-1)=nM(I, J-1)+nm4*M(I, J); % now do pushes to the
                nM(I+1, J-1)=nM(I+1, J-1)+nm2*M(I, J); % nM stack
                if i>0,
                    nM(I-1, J-1)=nM(I-1, J-1)+nm3*M(I, J);
                end;
                M(I, J)=0;
                I=I+1;
            end; % while loop on I and i
        end; % for loop on J and j
        M=nM; %next M starts with what was pushed onto it
    end % m loop
    Ku(n, 1:(n+1))=UT; Km(n, 1:(n+1))=MT;
    Ki(n, 1:(n+1))=IT; Kr(n, 1:(n+1))=RT; Kc(n, 1:(n+1))=CT;
end% n loop
Ku
Km