function retval = fcnradfluxR(y,t)
% usage: retval = fcnradfluxR(y,t)
% description: returns rhs for a radial parabolic equation for MOL
% The problem is a radially symmetric diffusion problem with right
% Dirichlet boundary condition, left Neuman. The left Neuman condition
% at r = r0 takes the form   -*pi*r0*D*N_r(r0,t) = q.   The PDE is
%  N_t = (1/r)*(D*r*N_r)_r + f(r,t,N), 0 <= r0 < r < R
% where f is specified below.  Note that the actual boundaries and values
% at r0 and R, have to be built into this function. Column vector y is 
% array of unknown node values of the solution N. The form of both y and
% retval is a column vector of the same size. It is assumed that the nodes
% are equally spaced from r0 to R and that the spacing size,
% dr, is such that dr < r0/2 (needed for validity of discretization).
% Note order of arguments: 
% for Matlab it has to be reversed from the octave/lsode order y,t.
% Also note: for this function to have accuracy at left endpoint,
% you must choose discretization such that dr < r0/2.

% test problem, r0>0: N(r,t)=exp(r^2/2-t)

% initialize boundaries and BCs
r0 = 1;
R = 2;
NR = exp(2-t);
% diffusion coefficient
D = 1;
% next, total amount of material crossing circle r=r0 into medium
q = -2*pi*D*r0*exp(0.5-t);
% end initializing
% flux term at left boundary 
Nrr0 = -q/(2*pi*r0*D);
% find spacing and rnodes assuming Dirichlet BCs
n = length(y);
rnodes = linspace(r0,R,n+1)';
dr = (R-r0)/n;
% redefine r0 to be at ficticious "known" left boundary
r0 = r0 - dr;
% redefine rnodes with the new r0
rnodes = [r0;rnodes];
% now approximate the value of N at new r0
Nr0 = y(2) - 2*Nrr0*dr; % this uses centered difference
rnodesh = linspace(r0+dr/2,R-dr/2,n+1)'; % half nodes of r
retval = [Nr0;y;NR];
% functional part of rhs...here (1+r*r)*N
fnval = -(3+rnodes.*rnodes).*retval;
% derivative part of rhs
retval(2:n+1) = D*(rnodesh(2:n+1).*(retval(3:n+2)-retval(2:n+1)) - ...
rnodesh(1:n).*(retval(2:n+1)-retval(1:n)))./(dr*dr*rnodes(2:n+1));
% add in functional part of rhs and trim
retval = retval + fnval;
retval = retval(2:n+1);