function [x,w]=sde1(afun,bfun,ts,x0)
% function [x,w]=sde1(afun,bfun,ts,x0)
% First-order accurate scheme to integrate the Ito SDE
%        dx = a(x,t)dt + b(x,t)dW
% for a realisation of one (scalar) Wiener noise W(t).
%
% x0   - column vector of initial condition
% afun - user supplied function afun(x,t) computes the
%        coefficient column vector of the dt term
% bfun - user supplied function bfun(x,t) computes the
%        coefficient column vector of the dW term
% ts   - row vector [t0 t1 ... tfin] of times at which
%        x(t) is computed using steps of size diff(ts)
% x    - columns are the solution at times given in ts
% w    - row of the Wiener process at times ts (w0=0)
%
% Ref: Kloeden & Platen, "Numerical solution of
% stochastic differential equations", Section 11.1 
%
% Corresponding Stratonovich SDE is dx=(a-bb'/2)dt+b.dW
%
% (c) Tony Roberts, 3/6/98, aroberts@usq.edu.au
%
d=length(x0);           % dimensionality of the problem
dt=diff(ts);            % vector of delta t
rdt=sqrt(dt);           % square root delta t
nt=length(dt);          % number of time steps
dw=randn(1,nt).*rdt;    % increments of the noise process
dd=(dw.^2-dt)./(2*rdt); % auxilary factors
%
x=zeros(d,nt+1);
x(:,1)=x0;
for n=1:nt
   a=feval(afun,x(:,n),ts(n));
   b=feval(bfun,x(:,n),ts(n));
   db=feval(bfun,x(:,n)+a*dt(n)+b*rdt(n),ts(n)) -b;
   x(:,n+1)=x(:,n)+a*dt(n)+b*dw(n)+db*dd(n);
end
%
w=cumsum([0 dw]);
end