What is stochastic differential equation and its need?
One way to approach your problem is to use an Extended Kalman Filter with a small step size.  $x,y,$ and $z$ are your state variables.  Runge Kutta works better than Euler.  The ODE solvers in Matlab also work, but need to be set up for a single small time step.  It looks like at each time step, you propagate the differential equations forward in time, and the perturb the new states with noise, and then propagate those, repeat, ...repeat....

Google ( Lorenz Equation Matlab ) for code examples.  You have a Noisy Lorenz system.  

$$ \left| \begin{array}{c} x[k+1] \\ 
y[k+1] \\ z[k +1] \end{array} \right|= \text{ODE STEP} \left(\left| \begin{array}{c} x[k] + \text{noise}_x \\ 
y[k] + \text{noise}_y \\ z[k]+ \text{noise}_z \end{array} \right|
\right)
$$

Noisey Van Der Pol Example,

This is an old homework problem with some unnecessary feature like measurement noise, but you can modify it to your needs.  

There is a way for the Matlab ODE solver to only iterate from 0 to h instead of calculating 41 sub intervals.  but I don't recall how after a few years

You can use a 4th order Runge Kukta 
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

there is also 
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_(SDE)

You should note that, given your problem statement, you don't really have a pure SDE, The Lorenz is a chaotic system, with a noise perturbation.

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.  

% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 1.2 $  $Date: 2002/06/17 13:20:38 $  

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];  

% HW 8  part c
clear all
h=.01;
tspan = [0:h:60];
y0 = [1; 0];
epsilon = [2 1 .2] ;
q=5;
r=1;  

y1=zeros(2,length(tspan));
F1=zeros(2,2,length(tspan));
y1(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(1));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y1(:,i-1));
    y1(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y2=zeros(2,length(tspan));
y2(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(2));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y2(:,i-1));
    y2(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y3=zeros(2,length(tspan));
y3(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(3));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y3(:,i-1));
    y3(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end  

figure(1)
plot(tspan,y1(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y1=y1(:,startp:endp);
t1=tspan(startp:endp);
h1=plot(y1(1,:),y1(2,:));
set(get(h1,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)')
title([' Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(1))])
%
%Plot of the second solution
 figure(2)
 plot(tspan,y2(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
 y2=y2(:,startp:endp);
t2=tspan(startp:endp);
 h2=plot(y2(1,:),y2(2,:));
set(get(h2,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(2))])
%
%  Plot of the third solution
figure(3)
 plot(tspan,y3(1,:))
[plot_limits]=ginput(2) %select 3  orbits
 startp=round(plot_limits(1,1)/h);
 endp=round(plot_limits(2,1)/h);
y3=y3(:,startp:endp);
 t3=tspan(startp:endp);
h3=plot(y3(1,:),y3(2,:));
set(get(h3,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['van der Pol Equation, \epsilon = ',num2str(epsilon(3))])
 z1=y1(1,:)+sqrt(r)*randn(size(y1(1,:)));  % add measurement noise
z2=y2(1,:)+sqrt(r)*randn(size(y2(1,:)));
z3=y3(1,:)+sqrt(r)*randn(size(y3(1,:)));