function [residue] = ce_res (Y,G,nlines,this_lambda)

this_size = length(this_lambda);
nresp  = size(G,1) / nlines;
nref   = size(G,2);

%  Compose V matrix
V = zeros(this_size,this_size);
for n = 1:this_size
  V(n,:) = (this_lambda.^(n-1))';
end

% Solve for residues
fprintf(1,'Solving for residues ...');
residue = zeros(nref*nresp,this_size);
for p  = 1:nref
  for k = 1:nresp
    % Solve for 2 residues for each DOF for each lambda
    H = Y(k:nresp:k+nresp*(this_size - 1),p);
    residue(k+(p-1)*nresp,:) = (V\H)';
  end
end

return