function [F4th,G4th,F5th,G5th] = rkf45_lin(A,B,h)
%
%  This function computes the F and G matrices of the 4th- and
%  5th-order approximations of Runge-Kutta-Fehlberg 4/5, when used
%  to simulate linear systems.
%
%         [F4th,G4th,F5th,G5th] = rkf45_lin(A,B,h)
%
%  Input parameters:
%  -----------------
%
%     A    :=  System matrix of continuous-time system.
%
%     B    :=  Input matrix of continuous-time system.
%
%     h    :=  Desired step length, can be either positive or negative.
%
%  Output parameters:
%  ------------------
%
%     F4th :=  System matrix of discrete-time system,
%              4th-order approximation.
%
%     G4th :=  Input matrix of discrete-time system,
%              4th-order approximation.
%
%     F5th :=  System matrix of discrete-time system,
%              5th-order approximation.
%
%     G5th :=  Input matrix of discrete-time system,
%              5th-order approximation.
%
I = eye(size(A));
A0 = A;
B0 = B;
%
F1 = I + (1/4)*h*A0;
G1 = (1/4)*h*B0;
A1 = A*F1;
B1 = A*G1 + B;
%
F2 = I + (3/32)*h*A0 + (9/32)*h*A1;
G2 = (3/32)*h*B0 + (9/32)*h*B1;
A2 = A*F2;
B2 = A*G2 + B;
%
F3 = I + (1932/2197)*h*A0 - (7200/2197)*h*A1 + (7296/2197)*h*A2;
G3 = (1932/2197)*h*B0 - (7200/2197)*h*B1 + (7296/2197)*h*B2;
A3 = A*F3;
B3 = A*G3 + B;
%
F4 = I + (439/216)*h*A0 - 8*h*A1 + (3680/513)*h*A2 - (845/4104)*h*A3;
G4 = (439/216)*h*B0 - 8*h*B1 + (3680/513)*h*B2 - (845/4104)*h*B3;
A4 = A*F4;
B4 = A*G4 + B;
%
F5 = I - (8/27)*h*A0 + 2*h*A1 - (3544/2565)*h*A2 + (1859/4104)*h*A3 - (11/40)*h*A4;
G5 = -(8/27)*h*B0 + 2*h*B1 - (3544/2565)*h*B2 + (1859/4104)*h*B3 - (11/40)*h*B4;
A5 = A*F5;
B5 = A*G5 + B;
%
F4th = I + (25/216)*h*A0 + (1408/2565)*h*A2 + (2197/4104)*h*A3 - (1/5)*h*A4;
G4th = (25/216)*h*B0 + (1408/2565)*h*B2 + (2197/4104)*h*B3 - (1/5)*h*B4;
F5th = I + (16/135)*h*A0 + (6656/12825)*h*A2 + (28561/56430)*h*A3 - (9/50)*h*A4 + (2/55)*h*A5;
G5th = (16/135)*h*B0 + (6656/12825)*h*B2 + (28561/56430)*h*B3 - (9/50)*h*B4 + (2/55)*h*B5;
%
return