Looks correct! No issue because I checked with the built-in @odephas2 function to generate the 2-D phase plane plot.
Edit: Upon rechecking your original code, I found that there were no errors during execution. However, I discovered that the code within the second_order_ode() function was incorrect. Please refer to the comment within the code for the necessary corrections.
ODE:
State-space:
function dydt = second_order_ode(t, y, param)
gamma - alpha*y(2) - beta*y(1)];
opts = odeset('OutputFcn',@odephas2, 'Stats', 'on');
[t, y] = ode45(@(t, y) second_order_ode(t, y, param), tspan, y0, opts);
22 successful steps
0 failed attempts
133 function evaluations