% kubo oscillator
clear, close all
Nmax = 1e5;
N = 1e3;
dt = 0.01;
Tf = N*dt;
x = zeros(N+1,1);
y = x;
tt = linspace(0,Tf, 1001);
dW = sqrt(dt)*randn(N,1);
x(1,1) = 0;
y(1,1) = 1;
for ii = 1:N
    x(ii+1,1) = x(ii,1) - 0.5*x(ii,1)*dt + y(ii,1)*dW(ii,1);
    y(ii+1,1) = y(ii,1) - 0.5*y(ii,1)*dt - x(ii,1)*dW(ii,1);
end
 plot(tt,x,tt,y)
    