function [x,w]=sde1m(abfun,ts,x0,pind)
% function [x,w]=sde1m(abfun,ts,x0,pind)
% First-order accurate scheme to integrate the Ito SDE
%        dx = a(x,t)dt + b(x,t)dW
% for a realisation of m>1 Wiener noises W(t).
%
% abfun - user supplied function [a,b]=abfun(x,t) 
%         computes the coefficient column vector a(x,t)
%         of the dt term, and the coefficient matrix b(x,t) 
%         of the dW term (m columns)
% ts   - row vector [t0 t1 ... tfin] of times at which
%        x(t) is computed using steps of size diff(ts)
% x0   - column vector of initial condition
% pind - controls approx of multiple stochastic integrals,
%        larger is better but more time consuming,
%        recommended minimum is max(diff(ts)), 0 is possible
% x    - columns are the solution at times given in ts
% w    - m rows of the Wiener processes at times ts (w0=0)
%
% Ref: Kloeden & Platen, "Numerical solution of
% stochastic differential equations", Section 11.1 
%
% (c) Tony Roberts, 3/6/98, aroberts@usq.edu.au
%
d=length(x0);           % dimension of state-space
[a,b]=feval(abfun,x0,ts(1));
m=size(b,2);            % number of independent noises
dt=diff(ts);            % vector of delta t
rdt=sqrt(dt);           % square root delta t
nt=length(dt);          % number of times
dw=randn(m,nt).*rdt(ones(m,1),:); % increments of the noises
%
x=zeros(d,nt+1);
x(:,1)=x0;
for n=1:nt
   %
   % compute the I matrix
   p=ceil(abs(pind)/dt(n));
   xi=dw(:,n)/rdt(n);
   mu=randn(m,1);
   if p>0,
      rr=1 ./(1:p);
      rrs=rr(ones(m,1),:);
      rhop=1/12-sum(rr.^2)/(2*pi^2);
      eta=randn(m,p);
      zet=randn(m,p);
      zetxi=(zet*rr')*xi';
      zetet=zet*(rrs.*eta)';
   else rhop=1/12; end
   muxi=(mu*xi')*sqrt(rhop);
   ip=xi*xi'/2;
   if p>=0, ip=ip +(muxi-muxi'); end
   if p>0, 
      ip=ip+( (zetxi-zetxi')/sqrt(2) +(zetet-zetet')/2 )/pi;
   end
   ip=dt(n)*(ip -diag(0.5*ones(m,1)) );
   %
   % supporting values 
   y=(x(:,n)+a*dt(n))*ones(1,m) +b*rdt(n);
   bb=zeros(d,1);
   for j1=1:m %sum 
      [ay,by]=feval(abfun,y(:,j1),ts(n));
	  bb=bb+(by-b)*ip(j1,:)';
   end
   %
   % time step
   x(:,n+1)=x(:,n)+a*dt(n)+b*dw(:,n) +bb/rdt(n);
   if n~=nt, [a,b]=feval(abfun,x(:,n+1),ts(n+1)); end
end
%
w=cumsum([zeros(m,1) dw]')';
end