You should collect your unknowns P and N in a vector x, and define a function for the right-hand side of your ODE system, something like this:
function dxdt = rhs(t,x,b,r_i,r_o,my_q)
P = x(1);
N = x(2);
dPdt = (b-r_i(N)-r_o(N))*P+r_i(N)*N;
dNdt = (b+my_q)*P-my_q*N;
dxdt = [dPdt;dNdt];
end
Save this as rhs.m.
Then create a script where you define all parameters and functions, as well as initial values. The ODE solvers do not take extra function parameters explicitly, so we define an inline function f of only t and x.
% Define parameters and functions going into your ODE system
b = 1;
r_i = @(N) sqrt(N) ;
r_o = @(N) exp(-N);
my_q = -4;
% Define an inline function f(t,x), where the extra parameters to rhs are
% "baked in":
f = @(t,x) rhs(t,x,b,r_i,r_o,my_q);
% Initialize:
P0 = 1;
N0 = 2;
x0 = [P0;N0];
% Solve:
[t,x] = ode45(f,x0,[0,100]);
P = x(:,1);
N = x(:,2);
plot(t,P);
If your system diverges very quickly, it may be useful to start with a very short time interval, to see what happens. For some parameter sets you system may be stiff, and ode45 may be very slow. If so, try e.g. ode15s.
