From the looks of your 3 coupled ordinary differential equations it seems unlikely that they have solutions that are neat elementary analytical functions. The easiest way to go about this problem then is to integrate them numerically. For that you should write one function that defines the three equations for dS/dt, dX/dt and dP/dt, and then use ode45 (or one of its siblings) to integrate that function from your initial conditions. See the help and documentation to ode45.
The ode-function should look something like this:
function dSdtdXdtdPdt = your3odes(t,SXP,other_pars)
D = other_pars(1);
Xi = other_pars(2);
Pi = other_pars(3);
% etc...
S = SXP(1);
X = SXP(2);
P = SXP(3);
dPdt = D*(Pi-P) + %the rest of that equation
dXdt = D*Xi - X* % the rest f that equation
dSdt = % same here
% Combine the 3 time-derivatives into one column-array
dSdtdXdtdPdt = [dSdt;dXdt;dPdt];
end
Then integrate that over your time-span of interest from your initial conditions [S0;X0;P0].
HTH
