diary
% script: radialNeumann.m
% the solution to this problem is N(r,t)=exp(r^2/2-t)
% get the initialization out of the way
n = 20;rnodes=linspace(1,2,n+1)';rnodes=rnodes(1:n);
% next the initial condition
y0 = exp(rnodes.*rnodes/2);
% calculate the solution
sln = lsode('fcnradfluxR',y0,[0,2]);
% plot solution and exact solution at t=2
clearplot;hold on;grid
plot(rnodes,sln(2,:)')
plot(rnodes,exp(rnodes.*rnodes/2-2))
% not as good as the Dirichlet problem, but this is typical
% check the maximum error
e1 = norm(sln(2,:)'-exp(rnodes.*rnodes/2-2),inf)
% ok, let's halve the stepsize
n = 40;rnodes=linspace(1,2,n+1)';rnodes=rnodes(1:n);
y0 = exp(rnodes.*rnodes/2);
sln = lsode('fcnradfluxR',y0,[0,2]);
plot(rnodes,sln(2,:)')
e2 = norm(sln(2,:)'-exp(rnodes.*rnodes/2-2),inf)
e1/e2
% once again, about 4.