function [t,X] = Platen2DSLVectorized(X0,A,g0,g,tspan,h,dW)
%Platen DSL scheme of weak order 2 to approximate the solution of 
% dX=[AX+g_0(X)]dt + g{1}(X)dW{1} on the time interval tspan with step size h,
% Wiener increments dW and initial value X0

%number of stochastic integrals
M=length(g);
if M>1
    error('The Platen method is only implemented for a single stochastic integral')
end

% Create temporal grid
t = tspan(1):h:tspan(2);
n = length(t);

%Initialise solution
X = repmat(X0,[1,size(dW{1},2)]);

% Apply exponential Platen 2.0

%Calculate exponentials
expDeltaL_n = expm(A*h); % e^{\Delta L_n}

%time stepping
for i = 2:n
    %Calculate stage values
    
    H1 = X;
    g0H1=g0(H1);
    g1H1=g{1}(H1);
    H2 = X + (g0H1*h + g1H1.*dW{1}(i-1,:));
    H3 = X + (g0H1*h + g1H1*sqrt(h));
    H4 = X + (g0H1*h - g1H1*sqrt(h));
    
    g1H3exp=g{1}(expDeltaL_n*H3);
    g1H4exp=g{1}(expDeltaL_n*H4);
    
    %Update solution
    X = expDeltaL_n*X + 0.5*(g0(expDeltaL_n*H2) + expDeltaL_n*g0H1)*h...
        + 1/4*(g1H3exp+ g1H4exp + 2*expDeltaL_n*g1H1).*dW{1}(i-1,:)...
        + 1/4*(g1H3exp - g1H4exp).* (dW{1}(i-1,:).^2 - h)/sqrt(h);
end