Dear Matlab experts or mathematicians:
I am a biology student and need to correctly code a diffusion equation (equation 6.43) shown in the snapshot below from the Mathematics of Diffusion by Crank. My problem is in how to correctly calculate and implement the roots and use them inside equation 6.43 below. Here is what I have done (assuming that I only need to know the first root of B cotB + L-1 =0):
First: the roots of equation 6.41 were calculated:
F = @(z) 1 - z.*cot(z) - L; % z is the root of the equation
z = linspace(0,3,1000);
G=F(z);
plot(z,F(z),'k','linewidth',2);
hold on;
plot(z,0*z,'k');
% by inspecting the plot a value (K) close to zero can be obtained
root= fzero(@(z) 1 - z.*cot(z) - L, K); % note that fzero fails at large L values ==> use fminbnd instead
end
now the first root is obtained which is by the way is ~ 3.14 since L > 1000. Next, the equation 6.43 is solved as follows:
D= 1e-18 %cm/s
a = 2e-7 %cm
tmax = 1E8; % in sec
t = linspace(0,tmax);
for n=1:10:1000
zeta=(n*root)^2; % I am not sure about if n*root is the correct implementation !!!
P=(zeta*(zeta + (L*(L-1))));
y=(6*L^2)*exp(-D*zeta*t/a^2)/P;
y=y+y;
w=1-y;
end
plot(t,w,'k','linewidth',2)
The result I am getting is not right and I guess that the way I have implemented the nth roots in equation 6.41 is incorrect. My maths background is not great and so I suspect that my problem is in understanding the mathematical formula. I mainly do not understand why I need to know the other roots if I know the first root solution for equation 6.41??
Could you please help me out.
Many thanks in advance.