function [F] = ff3(A,h,algor)
%
%  This function computes the F-matrix for different integration algorithms
%
%  1 = Nystrom-Milne 3rd order
%  2 = Nystrom-Milne 4th order
%  3 = Classical Milne algorithm
%
Ah = A*h;
Ah2 = Ah*Ah;
[n,n] = size(Ah);
Z = zeros(n);
I = eye(n);
%
%  1 = Nystrom-Milne 3rd order
%
if algor==1,
  F31 = Ah2/9;
  F32 = I + 2*Ah/3 - 2*Ah2/9;
  F33 = 4*Ah/3 + 7*Ah2/9;
  F = [ Z , I , Z ; Z , Z , I ; F31 , F32 , F33 ];
  return,
end
%
%  4 = Nystrom-Milne 4th order
%
if algor==2,
  F41 = -Ah2/9;
  F42 = 4*Ah2/9;
  F43 = I + 2*Ah/3 - 5*Ah2/9;
  F44 = 4*Ah/3 + 8*Ah2/9;
  F = [ Z, I, Z, Z; Z, Z, I, Z; Z, Z, Z, I ; F41, F42, F43, F44 ];
  return,
end
%
%  3 = Classical Milne algorithm
%
if algor==3,
  F41 = Ah/3;
  F42 = 8*Ah2/9;
  F43 = I + Ah/3 - 4*Ah2/9;
  F44 = 4*Ah/3 + 8*Ah2/9;
  F = [ Z, I, Z, Z; Z, Z, I, Z; Z, Z, Z, I ; F41, F42, F43, F44 ];
  return,
end
%
return
