%
clear all

N=1001;
Dt = 1;
t = Dt*[0:N-1]';

x = zeros(N,1);
y = zeros(N,1);
f = zeros(N,1);

Np = floor(N/10);
p = random('unid',N,Np,1);
cexp = 0.5;
x(p) = random('exp',cexp,Np,1);

figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 2 ] );
plot(t,x,'k-','LineWidth',2);

cf = 0.2;
for i=[2:N]
    f(i) = cf*(y(i-1)^2);
    dydt = x(i) - f(i);
    y(i) = y(i-1) + dydt * Dt;
    if( y(i) < 0 )
        y(i)=0;
    end
end

subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 1.1*max(f) ] );
plot(t,f,'k-','LineWidth',2);

subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 1.1*max(y) ] );
plot(t,y,'k-','LineWidth',2);

dlmwrite('nonlinearresponse.txt', [t,f,x], 'delimiter', '\t' );


