function [H0 H1] = solvemodel(A0,A1,A2,A3,A4,F)
% solve_model.m
% Uses approach set out in Simple_Example.tex to solve model
%
% A1*Z(t) = A0*Z(t-1) + A2*Z(t+1) + A3*F*z(t) + A4*v(t)
%
% Equilibrium is
%      Z(t) = H0*Z(t-1 + H1*v(t)
%
% where
%      H0 = (A1 - A2*H0 -A3*F)\A0 and
%
%      H1 = (A1 - A2*H0 -A3*F)\A4.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H0 = zeros(size(A1));
H1 = zeros(size(A1));

tolerance = 1e-7;
t=1;
displayevery = 1000;
maxiter = 10000;
crit(1) = 1;

disp(' ');disp('... computing the fixed point of [H0] ...');
while crit(t) > tolerance
    H0update = (A1 - A2*H0 -A3*F)\A0;    
    H0diff = abs(H0update - H0) ;
    H0 = H0update ; 
    t = t+1;
    crit(t) = max(norm(H0diff,1));
    
    if (rem(t,displayevery) == 0)
        outmess = ['t= ',int2str(t), '       crit: ', num2str(crit(t)), '    tolerance: ', num2str(tolerance)];
        disp(outmess);
    end%if (rem(t...
    if t > maxiter
        crit(t) = 0; % to force end of loop
        disp('max iteration reach before convergence')
    end
end

disp(['convergence achieved at t= ',int2str(t)]);disp(' ');

H1 =  (A1 - A2*H0 -A3*F)\A4; 