function G3_new = unfoldg3(G3,n)

% This function is designed to "unfold" matrix g_3 generated in the third
% order approximation of a policy function in Dynare.
% Arguments:
% G3 => matrix g_3, normally stored at oo_.dr.g_3, after the computation of
% a third-order approximation;
% n => number of state variables + shocks in the model

% Step 1: partition G3 matrix in "n" matrices with number of columns equal
% to n, n-1, n-2, n-3... 1.

i = 1;
j = 1;
nocols = 0;

while j < n+1
while i < n+1
    eval(['G3' num2str(i) num2str(j) '= G3(:, nocols + (- n + (n + 1.5)*i - (i^2)/2) :  nocols + (i*n - (0.5*(i^2) - (0.5*i))) );']);
    if j>i
        eval(['G3' num2str(i) num2str(j) '= [];']);
    end
    i = i+1;
end
for k = j+1:i-1
    nocols = nocols + size(eval(['G3' num2str(k) num2str(j)]),2);
end
i = 1;
j = j+1;
end

% Step 2: aggregate the partial matrices to form G31, G32, ... in such a way the problem can be solved using unfoldg2.m.

for i = 1:n
    eval(['G3' num2str(i) ' = [];']);
end

ncol1 = 0;
for i = 1:n
for j = 1:n
    eval(['G3' num2str(i) ' = [G3' num2str(i) ' G3' num2str(j) num2str(i) '];'])
    if i == 1
        ncol1 = ncol1 + size(eval(['G3' num2str(j) num2str(i)]),2);
    end
end
if size(eval(['G3' num2str(i)]),2) < ncol1
   subG = [];
   for n1 = 1:i-1
       for n2 = 1:n-(n1-1)
           eval(['subG = [subG G3' num2str(n1) '(:,i + n*(n1-1) +(n2-1)*n)];']);
       end
   end
   eval(['G3' num2str(i) '= [subG G3' num2str(i) '];']);
   eval(['G3' num2str(i) ' = unfoldg2(G3' num2str(i) ',n);']);
else
    eval(['G3' num2str(i) ' = unfoldg2(G3' num2str(i) ',n);']);
end
end

% Step 3: aggregate matrices G31, G32, ... to obtain G3_new

G3_new = [];
for i = 1:n
    G3_new = [G3_new eval(['G3' num2str(i)])];
end