I'm not usually in the habit of solving other people's homework, but today's your lucky day I guess.
So, you have a system of coupled ordinary differential equations:
α'(t) = s(β-α(β+1)-qα²)
β'(t) = (-β-αβ+γ)/s
γ'(t) = w(α-γ)
and you want to solve for
y = [α(t) β(t) γ(t)]
with 0 < t < 10, s = 1, q = 1, w = 0.1610. The way to do that in matlab is to define a function that computes the derivative ([α'(t) β'(t) γ'(t)]), and throw this in one of the ODE solvers (ode45
is a good first bet):
s = 1;
q = 1;
w = 0.1610;
% Define y(t) = [α(t) β(t) γ(t)] = [y(1) y(2) y(3)]:
deriv = @(t,y) [...
s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2) % α'(t)
(-y(2) - y(1)*y(2) + y(3))/s % β'(t)
w * (y(1)-y(3)) % γ'(t)
];
% initial value
y0 = [30 1 10];
% time span to integrate over
tspan = [0 10];
% solve ODE numerically
[t, y] = ode45(deriv, tspan, y0)
This will output
y =
30.0000 1.0000 10.0000
28.5635 0.9689 10.0049 % numerical solutions at times t
27.2558 0.9413 10.0094
26.0603 0.9166 10.0136
... ... ...
= α(t) = β(t) = γ(t)
t =
0
0.0016
0.0031 % corresponding times
0.0047
...
You can plot all this like so:
figure, clf, hold on
plot(t, y(:,1), 'r')
plot(t, y(:,2), 'g')
plot(t, y(:,3), 'b')
legend('\alpha(t)', '\beta(t)', '\gamma(t)')
which results in this figure: