Milstein's high order method graph stays at 0

I'm trying to graph the following SDE: dX = rX(K-X)dt +beta*sqrt(X)*dW using the following code:

rand('state',100);
beta = 0.25;
r = 1; K = 10; X0 = K; %problem parameters
T = 1; N = 2^11; dt = T/N;
dW = sqrt(dt)*randn(N,N);
Xmil = zeros(N,1);
for i = 1:N
    Winc = sum(dW(:,i));
    Xmil = Xmil + dt*r*Xmil.*(K-Xmil) + beta*(Xmil.^0.5).*Winc + 0.5*(beta^2)*(Xmil.^0.5).*(Winc.^2 - dt);
end
Dt = zeros(N,1); %building an array for the time values we have taken
for j = 2:N
   Dt(j) = j*dt; 
end
plot(Dt, Xmil);
title('Solution of The Equation using Milstein Method')
xlabel('t');ylabel('X(t)');

however the graph keeps staying at 0 for every value of beta I'm inputing.

[Amit R]

I've managed to find the solution on my own, needed to change up the initial values to non zero and fix up a few kinks in the code based on an article I had.

for those interested the final code is as follows:

rand('state',100)
r = 2; K = 1; Xzero = 0.5; % problem parameters
T = 1; N = 500; dt = T/N;
beta = 0.25;

dW = sqrt(dt)*randn(N,N); % Brownian increments
Xmil = zeros(N,1); % preallocate array
Xtemp = Xzero*ones(N,1);
for j = 1:N
Winc = sum(dW(:,j),2);
Xtemp = Xtemp + dt*r*Xtemp.*(K-Xtemp) + beta*Xtemp.*Winc + 0.5*beta^2*Xtemp.*(Winc.^2 - dt);
end
Xmil = Xtemp;
Dt = zeros(N,1); %building an array for the time values we have taken
for j = 2:N
   Dt(j) = j*dt; 
end
plot(Dt,Xmil)
title('Solution of The Milstein')
xlabel('t');ylabel('X(t)');

Answer 1

[Alan Stevens]

Your line

 Xmil = Xmil + dt*r*Xmil.*(K-Xmil) + .....

probably needs to be

 Xmil(i+1) = Xmil(i) + dt*r*Xmil(i).*(K-Xmil(i)) + .... etc.

However, if this is the case, you will still need a non-zero initial value for Xmil, or all terms will inevitably stay at zero!

[Amit R]

Thanks! I managed to fix this on my own but didn't update this post.