Summation of series with two index m,n

Hi,

I'm trying to simulate the following math equation( please see attached file name as equation).

There is term name is d1 which is my main concern it has two index in a single summation m and n

I want to verify that I have written right code or not

Please have a look on the attached file name as equation and in that, d1 is given which has a summation over m and in that summation it has a different index n

I put d1=0 after each calucaltion of phi for each n is it right?

In short d1 is given in equation file and I want to know that have I written the correct code of phi1 (as given in equation file) ?

Code is:

q=1.602e-19;
h=8.85e-12;
Eox=7.8*h;
Esi=11.8*h;
Ni=1.1e16;
Na=1e16;
%for Vgs=0:0.1:0.5
Vgs=0.1;
L=45e-9;
tox=3e-9;
tsi=10e-9;
Vds=1;
Q=-q*Na;
V=Vgs;
t1=0.02586;%thermal voltage in V%
Vbi=1.11/2;%source side potential%
phi=0;
d1=0;
phi_1=0;
%syms x y;
y=tox+tsi/2;
x=0:1e-10:L;
for n=1:1:20
    kn=n*pi/L;
    X1=((cosh(kn*(tox+tsi-y)))+(cosh(kn*(tox-y))))/(sinh(kn*tsi));
    %X_1=((cosh(kn*(tox+tsi-y_1)))+(cosh(kn*(tox-y_1))))/(sinh(kn*tsi));
    d2=L*(V*(1-(-1)^n))/((cosh(kn*tox))*pi*n);
    d4=1/(kn*Esi*tanh(kn*tsi))+1/(kn*Esi*sinh(kn*tsi))+tanh(kn*tox)/(kn*Eox) ;
    d3=Vbi*(1-(- 1)^n)/kn+(Vds*(-1)^(n+1))/kn+Q*(1-(-1)^n)/(Esi*kn*kn*kn);;
    for m=1:1:10
        k1=(2*m-1)*pi/(2*tox);
        A1=V/k1+(sin(k1*tox)*(Vbi-Vgs))/(tox*k1^2);
        A2=(V)/k1+(sin(k1*tox)*(Vbi+Vds-Vgs))/(tox*k1^2);
        d1=d1+(2*(sin(k1*(tox)))*kn*(A1+A2*((-1)^(n+1)))/(tox*(kn^2+k1^2)));
    end
    Dsf=((d1+d2-d3))/d4;
    phi=phi+2*((sin(kn*x))*X1*Dsf)/(kn*Esi*L);
    d1=0;
end
hold on
phi1=phi+Vbi+(x/L)*Vds+(Q*(x*L-(power(x,2))))/(2*Esi);
plot(x,phi1, 'b')