Implementing a diffusion equation (ions in sphere)...issue implementing nth roots

Dear Matlab experts:

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:

if true
      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:

if true

    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)
  end

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. Could you please help me out.

[Imam Gunapy]

In case equation 6.43 is asking for all the roots (rather than only the first root) of equation 6.41, then the problem becomes how to find all roots at a large L value and then to pass them from equation 6.41 to equation 6.43?

[David Goodmanson]

Hi Imam,

The equation does use more roots than the first one, but if L is large then your approximation of nth root = n*(first root) is not too bad. But there are a couple of other things that need to be done first. I believe that your expression for P is incorrect. And when you want to add up terms, the usual pattern is

y = 0;
for n = whatever
 yterm = (some calculation)
 y = y + yterm;
end

then in your case w = 1-y after that.

[Imam Gunapy]

Hi David, Thanks for your help. Could you please explain where the expression P is wrong. I am not able to spot the mistake.

[David Goodmanson]

You are using zeta in place of beta, and if you expand it out you get

 P = zeta^2 + zeta*L*(L-1);

which is not what the equation says.