Recently I had to accomplish the same task.
The "issue" is that ode45 of course changes the time step size due to the stability criteria. It implements a adaptive time step size.
This means that you need values from your stochastic process according to this adaptive time steps.
An easy way is to interpolate between the discretized values of your stochastic process and the demanded time of the ode45 function.
For demonstration I used an example from the Matlab ode45 documention, and applied a brownian motion to this model that can be seen as an artificial parameter.
So to combine those two models, this is what I've done:
tspan = [0 5];
y0 = 0;
W = [0; cumsum(randn(length(tt)-1,1))]/sqrt(length(tt)-1);
W = W*sqrt(tt(end));
B = b*tt' + sigma*W;
title([int2str(length(tt)) '-step version of Brownian motion and its mean'])
xlabel(['Drift ' num2str(b) ', diffusion coefficient ' num2str(sigma)])
[t,y] = ode45(@(t,y) vdp1(t,y,tt,B),tspan,[2; 0]);
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
Now this is the crucial part, for every call of ode45 to the function vdp1, a value from the stochastic process needs to be interpolated.
function dydt = vdp1(t,y, tt, B)
intB = interp1(tt,B,t);
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)-intB];
This is the solution of the ordinary 2nd order van der Pol equation (seen in the matlab documentation example):
And this is the solution where a stochastic process has induced some changes in every time step:
Copyright 1984-2014 The MathWorks, Inc.