More like this perhaps, though the end result is rather boring as your forcing function, F, is all zeros and your initial conditions are all zero!
% A*x' + B*x = F
% x' = A\(-B*x+F)
% Initial Conditions
h = 0;
a = 0;
h_dot = 0;
a_dot = 0;
x_vector = [h; a; h_dot; a_dot];
time = [0 5];
[T,X] = ode45(@fun, time, x_vector);
plot(T, X(:,1));
hold on
plot(T, X(:,2));
hold off
function dXdt = fun(~, x)
m = 1;
x_m = 0.05;
e = 0.4;
zeta = 0.1;
w_h = 0.5;
w_a = 1.0;
I_a = 0.25;
D_a = 0.2;
U = 5;
a = x(2); h_dot = x(3);
L = D_a*U^2*(a+h_dot/U);
M = [1 x_m; (m*x_m)/I_a 1];
C = [2*zeta*w_h 0; 0 2*zeta*w_a];
K = [w_h^2 0; 0 w_a^2];
G1 = [-L/m; (L*e)/I_a];
G2 = [0;0];
I = [1 0; 0 1];
Z = [0 0; 0 0];
A = [0.5*C M; I Z];
B = [K 0.5*C; Z -I];
F = [G1; G2];
dXdt = A\(-B*x + F);
end
