function [t,X] = Platen15Vectorized(X0,g0,g,tspan,h,dW,dZ)
%Platen scheme of strong order 1.5 to approximate the solution of 
% dX=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,2)]);

%do time stepping
for i = 2:n
    %Calculate stage values
    % H1
    H1 = X;
    g0H1 = g0(H1);      g1H1 = g{1}(H1);
    % H2
    H2 = X + g0H1*h + g1H1*sqrt(h);
    g0H2 = g0(H2);      g1H2 = g{1}(H2);
    % H3
    H3 = X + g0H1*h - g1H1*sqrt(h);
    g0H3 = g0(H3);      g1H3 = g{1}(H3);
    % H4
    H4 = H2 + g1H2*sqrt(h);
                        g1H4 = g{1}(H4);
    % H5
    H5 = H2 - g1H2*sqrt(h);
                        g1H5 = g{1}(H5);
    % update solution
    
    X = X + g1H1 .* dW{1}(i-1,:)...
                + 1/2 * (...
                            g0H2 - g0H3...
                                ).* (dZ(i-1,:) / sqrt(h))...
                + 1/4 * (...
                            g0H2 + 2*g0H1 + g0H3...
                                ) * h...
                + 1/4 * (...
                            g1H2 - g1H3...
                                ).* (dW{1}(i-1,:).^2-h)/sqrt(h)...
                + 1/2 * (...
                            g1H2 - 2*g1H1 + g1H3...
                                ) .* (dW{1}(i-1,:)*h - dZ(i-1,:))/h...
                + 1/4 * (...
                            g1H3 - g1H2 + g1H4 - g1H5...
                                ) .* (1/3 * dW{1}(i-1,:).^2 - h) .* dW{1}(i-1,:) / h;
end