function out=findboundary(model,data,plotoption)

l = length(model);   %number of items
binnum = ceil(l*.1); %select the number of bins to be 10% of the total number of points

diff = (data-model)./model; %(o-e)/e

[n xout] = hist( diff,binnum);  

sdd = std(diff);
maxn = max(n);
meanx = mean(xout);
temp1 = ones(maxn)*sdd*2;
temp2 = ones(maxn)*sdd*3;
%sdd = sdd+meanx;


if plotoption == 'p'
    subplot(2,2,1);
    hist( diff,binnum);  %plot histogram
    hold on;
    plot(temp1,1:maxn,'Color','black');
    plot(temp2,1:maxn,'Color','red');
    legend(['\mu_x=',num2str(meanx), ', \sigma_{diff}=',num2str(sdd)],'2\sigma_{diff}','3\sigma_{diff}');
    hold off;
    
    subplot(2,2,2);
    normplot(diff);

    subplot(2,2,[3,4]);
    plot(data-model);
    xlabel('day number');
    ylabel('difference data-model');
end;

ind_0 = find(n == 0);
%ind_pos = min(find(ind_0>0));
%if isempty(ind_pos)
%    out = sdd*3;
%else
%    out = xout(ind_0(ind_pos));
%end;
out = sdd*3;


