function [F] = ff(A,h,algor)
%
%  This function computes the F-matrix for different integration algorithms
%
%  1 = FE-BE predictor-corrector algorithm with one corrector
%  2 = FE-BE predictor-corrector algorithm with two correctors
%  3 = FE-BE predictor-corrector algorithm with three correctors
%  4 = FE-BE predictor-corrector algorithm with four correctors
%  5 = FE-BE predictor-corrector algorithm with five correctors
%
I = eye(size(A));
Ah = A*h;
Ah2 = Ah*Ah;
Ah3 = Ah2*Ah;
Ah4 = Ah3*Ah;
Ah5 = Ah4*Ah;
Ah6 = Ah5*Ah;
%
%  1 = FE-BE predictor-corrector algorithm with one corrector
%
if algor==1,
  F = I + Ah + Ah2;
  return,
end
%
%  2 = FE-BE predictor-corrector algorithm with two correctors
%
if algor==2,
  F = I + Ah + Ah2 + Ah3;
  return,
end
%
%  3 = FE-BE predictor-corrector algorithm with three correctors
%
if algor==3,
  F = I + Ah + Ah2 + Ah3 + Ah4;
  return,
end
%
%  4 = FE-BE predictor-corrector algorithm with four correctors
%
if algor==4,
  F = I + Ah + Ah2 + Ah3 + Ah4 + Ah5;
  return,
end
%
%  5 = FE-BE predictor-corrector algorithm with five correctors
%
if algor==5,
  F = I + Ah + Ah2 + Ah3 + Ah4 + Ah5 + Ah6;
  return,
end
%
return
