How can I write an SDE in Matlab?

You absolutely cannot use ode45 to simulate this. ODEs and SDEs are different beasts. Also note that rand in Matlab is uniformly, not normally, distributed, so you're also not even generating proper Wiener increments. Before proceeding, you should take some time to read up on SDEs and how to solve them numerically. I recommend this paper, which includes many Matlab examples:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

The URL to the Matlab files in the paper won't work; use this one.

If you are familiar with ode45 you might look at my SDETools Matlab toolbox on GitHub. It was designed to be fast and has an interface that works very similarly to Matlab's ODE suite. Here is how you might code up your example using the Euler-Maruyma solver and anonymous functions:

a = ...
b = ...
c = ...
f = @(t,x)[ x(1)*(x(2)-a)+x(3);
            1-b*x(2)-x(2)^2;
           -x(1)-c*x(3)]; % drift function
g = 10; % or @(t,x)10;    % diffusion function, constant in this case
dt = 1e-3;                % time step
t = 0:dt:1;               % time vector
x0 = [...];               % 3-by-1 initial condition
opts = sdeset('RandSeed',1);  % Set random seed
y = sde_euler(f,g,t,x0,opts); % Integrate
figure;
plot(t,y);

In your example, the drift function, $g(t,x)$, doesn't depend on the state variable, $x$. This is so-called additive noise. In this case, more complex SDE integration schemes like Milstein reduce to Euler-Maruyma. More importantly, if you do use more complex diffusion functions (e.g., multiplicative noise where $g(t,x) = g_1(t)x$ like in the equation at the beginning of your question), you need to specify what SDE formulation you're using, Itô (most common in finance) or Stratonovich. My SDETools currently defaults to Stratonovich, but you can change it via opts = sdeset('SDEType','Ito');.


Unfortunately, I don't think you can use ode45 to solve it, whereas, since $dW$ distributes like $N(0,\sqrt{dt})$ you must use the following discretisation: $$ x_{n+1,i}=x_{n,i}+f_i(x_{n,i})\Delta t+g_i(x_{n,i})x_{n,i}\sqrt{\Delta t}\cdot\epsilon_i, $$ $$ \epsilon_i=\mbox{randn()} $$ $$ x_{n,i}\approx x_i(t) $$