function out=handToyModel(data,effect)
steps=1000;
dt=0.01;
stopThreshold=0.001;
numberReps=1000;
stimulusRange=0.8:0.01:1.2;
noise=0.01;
numberStimuli=length(stimulusRange);

% effect=;

resultHist=zeros(numberStimuli,numberReps);

h=waitbar(0);

for j=1:numberStimuli
waitbar(j/numberStimuli,h);
    for i=1:numberReps
        
        L=0;
        R=0;
        LHist=[]; LHist(1)=L;
        RHist=[]; RHist(1)=R;
        
        Lin=1*effect+randn*noise;
        Rin=stimulusRange(j)*effect+randn*noise;
        
        Lin=Lin;
        Rin=Rin;
        
        LinHist=[]; LinHist(1)=Lin;
        RinHist=[]; RinHist(1)=Rin;
        
        alpha=0.01;
        beta=0.02;
        gamma1=0.01;
        bias=data;
        gamma2=gamma1*bias;
        
        gamma1=gamma1;
        gamma2=gamma2;
        
        delta=0.03;
        
        status=1;
        t=2;
        
        while status==1
            
            dLin=-delta*LinHist(t-1);
            dRin=-delta*RinHist(t-1);
            dL=-delta*LHist(t-1) + alpha*LinHist(t-1) + beta*LHist(t-1) - gamma1*RHist(t-1);
            dR=-delta*RHist(t-1) + alpha*RinHist(t-1) + beta*RHist(t-1) - gamma2*LHist(t-1);
            
            LinHist(t)=LinHist(t-1)+dLin;
            RinHist(t)=RinHist(t-1)+dRin;
            LHist(t)=LHist(t-1)+dL;
            RHist(t)=RHist(t-1)+dR;
            
            condition1=abs(LHist(t)-LHist(t-1))<stopThreshold;
            condition2=abs(RHist(t)-RHist(t-1))<stopThreshold;
            condition3=RHist(t)*LHist(t)<0;
            if  (condition1 & condition2 & condition3) | t>1000
                status=0;
                resultHist(j,i)=RHist(t)-LHist(t);
            else
                t=t+1;
            end
        end
    end
end
close(h);
out = resultHist;

% figure;hold on;
% plot(LinHist,'c');
% plot(RinHist,'m');
% plot(LHist,'b');
% plot(RHist,'r');