function diffusion_gradm

%% simulate a simple morphogen diffusion model in a ring of cells

%% model parameters

n=50; %number of cells

diffusion=5; %diffusion coefficient
deg=0.05; %degradation rate
source = 100; %source strength
thresh = 10;


%% simulation settings

tstart=0;  %start time
tend=200;   %end time
nrep=100;  %number of solution reports

%% initial conditions

init=zeros(n,1);

%% run the simulation -- do not change

diffop=circshift(eye(n),1)+circshift(eye(n),-1)-2*eye(n);
sol=ode45(@model_odes,[tstart,tend],init);
tout=linspace(tstart,tend,nrep);
stout=deval(sol,tout);
ms=max(max(stout))

%% plot the solution at nrep regular time intervals

for j=1:nrep
    ringplotter(n,stout(:,j),tout(j))
    pause(2/nrep)
end

%% space-time plots of solution

figure(21), clf
subplot(2,1,1)
imagesc(stout'), axis xy
colormap hot
colorbar
set(gca,'FontSize',14)
xlabel('cell number')
ylabel('time')
nt=get(gca,'YTickLabel')
set(gca,'YTickLabel',linspace(tstart,tend,length(nt)))

subplot(2,1,2)
bar(stout(:,end),'y')
axis([0 n+1 0 1.05*max(stout(:,end))])
set(gca,'FontSize',14)
xlabel('cell number')
ylabel('concentration')

%% specification of the ODEs

    function dydt=model_odes(t,y)
        dydt=[source;zeros(n-1,1)] - deg*y./(thresh + y) + diffusion*diffop*y;
    end

%% function for plotting solutions on a ring of cells -- do not change

function ringplotter(n,d,t)

%given data d on a ring of n cells
%plot value of d in each cell

a=2*pi/n;
r1=1;
r2=1.2;
figure(20)
clf
for i=1:n
    patch([r1*cos((i-1)*a),r2*cos((i-1)*a),r2*cos(i*a),r1*cos(i*a)],[r1*sin((i-1)*a),r2*sin((i-1)*a),r2*sin(i*a),r1*sin(i*a)],(d(i)/ms)*[1,1,1])
end
colormap gray, colorbar('SouthOutside')
axis tight, axis off
set(gca,'FontSize',16)
set(gca,'DataAspectRatio',[1 1 1])
title(['Morphogen concentration. Time = ',num2str(t)])

end
end