function retval = fcnradrR(y,t)
% usage: retval = fcnradrR(y,t)
% description: returns rhs for a radial parabolic equation for MOL
% The problem is a radially symmetric diffusion problem with 
% Dirichlet boundary conditions. The PDE of the problem 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 in this function. Column vector y is array
% of interior 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, both known. Note order of arguments: for
% Matlab it has to be reversed from the octave/lsode order y,t.

% test problem, r0>0: N(r,t)=exp(r^2/2-t)
% another test, r0=0: N(r,t)=1+exp(-t)*cos(r)

% initialize boundaries and BCs
r0 = 1;
R = 2;
Nr0 = exp(1/2-t);
NR = exp(2-t);
% diffusion coefficient
D = 1;
% end initializing
% find spacing and rnodes assuming Dirichlet BCs
n = length(y);
rnodes = linspace(r0,R,n+2)';
dr = (R-r0)/(n+1);
rnodesh = linspace(r0+dr/2,R-dr/2,n+1)'; % half nodes of r
retval = [Nr0;y;NR];
% functional part of rhs...here -(3+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);