main.m

%%% parameters
% forget what you know
clear
% set max N
N=500;
% set "turnover n", aka n0
n0=10;
% plotting parameters:
zmax=.2;
xmax=60;
% how many R values to use?
NRs=20;

%%% rate parameters
% set max decay
Rmax=1e-1;
Rmin=1e-2;
% set vector of decay rates
Rs=linspace(Rmin,Rmax,NRs);
% set leakiness
a0=0
a0=.1;

%%% creation rate (function)
% will loop over "n"
n=0:N;
% use a0 and n to set creation rate function
wp=a0+n.^2./(n0.^2+n.^2);



% LOOP!
for idx=1:length(Rs);
  % set decay function=R*n
  wm=Rs(idx)*n;
  % calculate steady state q
  q(idx,:)=qcalc_smart(wp,wm);
  % calculate <n>
  nbar(idx)=dot(n,q(idx,:));
end

%%% plot some stuff

if 0
% plot creation rate (NB: could do this BEFORE the loop but...)
figure(4)
plot(wp);
title('wp')
end

% plot nbar vs R
figure(1)
plot(Rs,nbar)
title('nbar vs Rs')

% plot q(R,n)
figure(2)
imagesc(q)
xlim([0 xmax])
title('q')

% replot q(R,n)
figure(3)
mesh(q)
title('q')
zlim([0 zmax])
xlim([0 xmax])
rotate3d on 
zlabel('p(n)')
ylabel('R index')
xlabel('n')

% plot q
figure(5)
plot(q','.-')
title('all the p')
xlim([0 xmax])
xlabel('n')
ylim([0 zmax])


% plot q for each R (a loop, with a pause)
figure(6)
for idx=1:length(Rs)
  plot(q(idx,:),'.-')
  title(sprintf('q, R=%g',Rs(idx)));
  xlim([0 xmax])
  pause(1)
end