usol =

As usual, we begin by identifying the critical point(s) of the system, if they exist. However, I believe your code is correct, although it may be less accurate due to the use of the Euler method. If you examine my results generated by the ode45 solver (with ), you will notice that the state rapidly increases to a very high value before eventually converging to the calculated critical point. It's also worth noting that I selected different initial values to avoid encountering a division-by-zero singularity.

syms u v

%% parameters

P = 0.75; % Production rate of u

b = 0.9; % Production rate of v stimulated by u

mu = 0.5; % Decay rate of u

k = 1; % Decay rate of v (constant in this scenario)

a = 9.9; % Different rates of a to test

%% find critical points

eq1 = P - a*u^2/v - mu*u == 0;

eq2 = b*u^2 - k*v == 0;

eqn = [eq1, eq2];

sol = solve(eqn);

ucrit = double(sol.u)

vcrit = double(sol.v)

%% call ode45 solver and plot result

tspan = [0 100];

x0 = [1; 0.5];

options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);

[t, x] = ode45(@(t, x) ode(t, x, P, b, mu, k, a), tspan, x0, options);

plot(t, x), grid on, legend('u', 'v'), xlabel('t')

%% Check if the states (u & v) settle near the critical point (equilibrium)

x(end,:)

%% Gierer-Meinhardt activator-inhibitor model

function dx = ode(t, x, P, b, mu, k, a)

u = x(1);

v = x(2);

dx(1,1) = P - a*u^2/v - mu*u;

dx(2,1) = b*u^2 - k*v;

end