% script: radialDirichlet.m
% the solution to this problem is N(r,t)=exp(r^2/2-t)
% get the initializations out of the way
n = 20;rnodes = linspace(1,2,n+1)';rnodes = rnodes(2:n); 
% next the initial condition
y0 = exp(rnodes.*rnodes/2);
 % calculate the solution
sln = lsode('fcnradrR',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))
% looks nice...let's see the max error
e1 = norm(sln(2,:)'-exp(rnodes.*rnodes/2-2),inf)
% not bad...let's test the quality of our solution
% since spatial error is O(dr^2), halving the stepsize
% should cut the error down by about 4. let's see..
n = 40;rnodes = linspace(1,2,n+1)';rnodes = rnodes(2:n); 
y0 = exp(rnodes.*rnodes/2);
sln = lsode('fcnradrR',y0,[0,2]);
% we won't bother with plots..
e2 = norm(sln(2,:)'-exp(rnodes.*rnodes/2-2),inf)
e1/e2
% well, it gives a reduction factor of about 4.