% code to simulate diffusion on a 2-d circular domain via Method of lines (FRAP)
function Diff_sim
 global N rj dr
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

 N = 50; % number of grid points 
 dr = 1/N; %spatial step size
 r = [0:N]/N; % grid points
 rj = r(2:end-1)'; % interior grid points
 
 u0 = ones(N-1,1); %initial data
  
tstep = 0.01; 
t_end = 1; %total time to run simulation
 
%specify the output points
tspan = [0:tstep:t_end];

 
% integrate the system of ode's:
[T,S] = ode23s(@deRHS,tspan, u0, odeset('maxstep',1));  

% plot the solution at each time step
  
figure(1)
plot(T,1-S(:,1), 'linewidth',2)
xlabel('t','fontsize',20)
ylabel('u(0,t)','fontsize',20)
if(cmykflg==1)
print('../../figs_c/chapt_5/FRAP_1','-deps','-cmyk')
end
 
 figure(2)
 plot(T,log(S(:,1)), 'linewidth',2)
xlabel('t','fontsize',20)
ylabel('ln(1-u(0,t))','fontsize',20)
if(cmykflg==1)
print('../../figs_c/chapt_5/ln_FRAP','-deps','-cmyk')
 end
p=polyfit(T(9:end),log(S(9:end,1)),1)

 
figure(3)
for j = 1:length(T)
    plot(r,[1-S(j,1),1-S(j,:),1])
    hold on
end
hold off

jlist = [1,8,15];
%make a contour color plot 
for jj = 1:3
    j = jlist(jj);
    tj=T(jj)
x=[-2:.05:2];
y=[-2:.05:2];

for kx = 1:length(x)
 for   ky=1:length(y)
    rsq = sqrt(x(kx)^2+y(ky)^2);
    
Sr(kx,ky) = evals(rsq,[r,10],[ S(j,1), S(j,:),0,0]);
end
end

[xx,yy] = meshgrid(x,y);
bottom = 0;
top  = 1;

figure(3+jj)
 
contourf(xx,yy,Sr,25 )
shading interp;
 caxis manual
caxis([bottom top]);
 end
 
%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global N rj dr
 
u = [s(1);s;0];
 up1= u(3:end);
um1=u(1:end-2);
 % here is the finite difference for the radially symmetric diffusion
 % operator
Fu = (((up1-s).*(rj+dr/2)-(s-um1).*(rj-dr/2))./rj)/dr^2;
 s_prime = Fu;

function sr = evals(r,R,S);
        % use linear interpolation to evaluate S at r
        jr = min(find(r<R));
         
        sr = S(jr-1)+(S(jr)-S(jr-1)).*(r-R(jr-1))./(R(jr)-R(jr-1));
        
