clear all;
% tau*dx/dt + x = D(t) -> dx = (D(t)*dt - x*dt)/tau
% initial condition: x(t=0) = x0;

% set up discritized time vector
dt = 1;
t = [0:dt:10];

% parameters of equation
a = 10;
D = a*ones(size(t));
tau = 2;
x0 = 0;

x(1) = x0; % initial condition
for i=2:length(t)
    dx = (D(i)*dt - x(i-1)*dt)/tau;
    x(i) = x(i-1) + dx;
end

% plot
plot(t,x,'r')
xlabel('time');ylabel('x');

% analytic solution
%at t=inf, dx/dt = 0 -> x=xinf = a
% x = xinf + (x0-xinf)e^(-t/tau);
xinf = a;
x_analytic = xinf + (x0-xinf)*exp(-t/tau);
hold on;
plot(t,x_analytic,'b')
legend('numeric','analytic')