clear all;clc;
nsteps = 500;
soln = zeros(2,nsteps+1); %allocate space for soln vector
soln(:,1)  = [1,2]; % whatever x0 and y0 initial conditions are

t(1) = 0;
for n = 1:nsteps
    zn = soln(:,n);   % solution from previous timestep
    D = 10;           % whatever parameter D is
    dt = 0.1;         % timestep size

    % nonlinear equation solver to find next time step value
    soln(:,n+1) = fsolve(@(znp1) (znp1(2)-zn(2))/dt + zn(2)/D - (znp1(1)-zn(1))/dt, soln(:,n) );
    t(n+1) = t(n) + dt; % increment time for a counter
end
plot(t,soln(1,:),t,soln(2,:)); legend('x(t)','y(t)')