Finite Difference Approach to Solving Schrödinger Equation

Particle in a box

n = ii = jj = 100;
l = 1;
Δ = l/n;
x0 = 0.;
x[i_]:= x0 + i Δ

h = Table[If[i==j,2./Δ^2,0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

vals[[n-1]]/π^2
vals[[n-2]]/π^2
vals[[n-3]]/π^2

Table[{i,vals[[n-i]]/(i π)^2},{i,1,n-1}]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

x[1]

vecs[[n-1]][[1]]
vecs[[1]][[1]]

num = 2
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}]]

Particle in a box (4th order algorithm)

n = ii = jj = 10;
l = 1;
Δ = l/n;
x0 = 0.;
x[i_]:= x0 + i Δ

h = Table[If[i==j,30./(12Δ^2),0]+ If[j==i+1 || j==i-1,-16/(12Δ^2),0]+If[j==i+2|| j ==i-2,1/(12 Δ^2),0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

vals[[n-1]]/π^2
vals[[n-2]]/π^2
vals[[n-3]]/π^2

Table[{i,vals[[n-i]]/(i π)^2},{i,1,100}]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

x[1]

vecs[[n-1]][[1]]
vecs[[1]][[1]]

num = 40
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}]]

Particle in a “complicated” box I

n = ii = jj = 1000;
l = 1;
Δ = l/n;
x0 = 0.;
x[i_]:= x0 + i Δ
u[x_]:= If[x<0.5,0,49 π^2]
Plot[u[x],{x,0,1}]

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

vals[[n-1]]/π^2
vals[[n-2]]/π^2
vals[[n-3]]/π^2

Table[{i,vals[[n-i]]/(i π)^2},{i,1,500}];

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

num = 3
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}]]

ListPlot[vecs[[n-11]]]

ListPlot[vecs[[n-2]]]

Particle in a triangular “box”

n = ii = jj = 1000;
l = 10;
Δ = l/n;
x0 = 0.;
x[i_]:= x0 + i Δ
u0 = 20;
u[x_]:= u0 x
Plot[u[x],{x,0,1}]

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

num = 1;
vals[[n-num]]
vals[[n-num]]/u0
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 5;
vals[[n-num]]
vals[[n-num]]/u0
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 15;
vals[[n-num]]
vals[[n-num]]/u0
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

Particle in a finite square well

n = ii = jj = 1000;
l = 1;
x0 = -1.;
xr = 1.;
Δ = (xr - x0)/n;
x[i_]:= x0 + i Δ
u0 = 51 π^2;
u[x_]:= If[-0.5<x<0.5,-u0,0]

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

Technique will only work for bound states, so find them.

boundVals = {};
boundVecs = {};
For[i=1,i≤ Length[vals],++i,If[vals[[i]]<0,{boundVals =Append[boundVals ,vals[[i]]],
boundVecs = Append[boundVecs,vecs[[i]]]}]]

(u0+boundVals)/π^2

Table[{i,vals[[n-i]]/(i π)^2},{i,1,500}];

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

num = 1;
(u0+boundVals[[num]])/π^2
ListPlot[Table[{x[i],boundVecs[[num]][[i]]},{i,1,Length[vecs[[num]]]}]]

num = 7;
(u0+boundVals[[num]])/π^2
ListPlot[Table[{x[i],boundVecs[[num]][[i]]},{i,1,Length[vecs[[num]]]}]]

NOTE : The range I chose (from -2 to +2) isn’t really big enough for the n=8 state illustrated below.

num = 8;
(u0+boundVals[[num]])/π^2
ListPlot[Table[{x[i],boundVecs[[num]][[i]]},{i,1,Length[vecs[[num]]]}]]

Harmonic oscillator

n = ii = jj = 1000;
xl = -20.;
xr = 20.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
u0 = 4.;
u[x_]:= x^2

x[n]

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

vals[[n-1]]

vals[[n-2]]

vals[[n-3]]

num = 1;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 10;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 150;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

Table[{i,vals[[n-i]]},{i,1,500}];

ListPlot[%,PlotRange->{{0,100},{0,400}}];
Plot[1 + 2x,{x,0,100}];
Show[%,%%]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

num = 150;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->{{-1,0},All}]

Compare minimum λ and Δ :

2 π/Sqrt[80.]

Harmonic oscillator (4th order)

n = ii = jj = 1000;
xl = -20.;
xr = 20.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
u0 = 4.;
u[x_]:= x^2

x[n]

h = Table[If[i==j,30./(12Δ^2)+u[x[i]],0]+ If[j==i+1 || j==i-1,-16/(12Δ^2),0]+If[j==i+2|| j ==i-2,1/(12 Δ^2),0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

vals[[n-1]]

vals[[n-2]]

vals[[n-3]]

num = 1;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 10;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

num = 150;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->All]

Table[{i,vals[[n-i]]},{i,1,500}];

ListPlot[%,PlotRange->{{0,100},{0,400}}];
Plot[1 + 2x,{x,0,100}];
Show[%,%%]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

num = 150;
vals[[n-num]]
Sqrt[vals[[n-num]]/u0]
ListPlot[Table[{x[i],vecs[[n-num]][[i]]},{i,1,Length[vecs[[n-num]]]}],PlotRange->{{-1,0},All}]

Compare minimum λ and Δ :

2 π/Sqrt[80.]

Inverted Gaussian

n = ii = jj = 2000;
xl = -10.;
xr = 10.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
u0 = 10000.;
u[x_]:= -u0 Exp[-x^2]

x[n]

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

{vals,vecs}=Eigensystem[h];

ListPlot[vals]

Technique will only work for bound states, so find them.

boundVals = {};
boundVecs = {};
For[i=1,i≤ Length[vals],++i,If[vals[[i]]<0,{boundVals =Append[boundVals ,vals[[i]]],
boundVecs = Append[boundVecs,vecs[[i]]]}]]

Breaks down when λ is on order of Δ

λ = 2/300.

Δ//N

boundVals

boundVecs[[1]];

Length[boundVals]

num = 3;
nn = Length[boundVals];
boundVals[[num]]
ListPlot[Table[{x[i],boundVecs[[num]][[i]]},{i,1,Length[boundVecs[[nn-num]]]}],PlotRange->{{-2,2},All}]

num = 82;
nn = Length[boundVals];
boundVals[[nn-num]]
ListPlot[Table[{x[i],boundVecs[[num]][[i]]},{i,1,Length[boundVecs[[num]]]}],PlotRange->{{-.5,.5},All},Joined->True]

Plot[u[x],{x,-10,10},PlotRange->All]

Compare minimum λ and Δ :

Δ
2 π/Sqrt[u0]

Particle in 2-D box

n = ii = jj = 51;
l = 1;
Δ = l/n;
x0 = 0.;
y0 = 0.;
x[i_]:= x0 + i Δ
y[i_]:= y0 + i Δ
ijToIndex[i_,j_]:= (i-1)*n +j    (* transform x-y coordinates to 1-D indexing *)
indexToij[i_]:= {Ceiling[i/(n-1)],i-(Ceiling[i/(n-1)]-1)*(n-1)} (* transform 1-D indexing to x-y *)
u[{x_,y_}]:= 0

NOTE : Tricky part is getting points on boundary right.

h = Table[If[i==j,4./Δ^2+ u[indexToij[i]],0]
              + If[j==i+1 && Mod[j,n-1]-1 ≠ 0,-1/Δ^2,0]
              + If[j==i-1 && Mod[i,n-1]-1 ≠ 0,-1/Δ^2,0]
              +If[(j== i+(n-1)||j== i-(n-1)),-1/Δ^2,0],{i,1,(n-1)^2},{j,1,(n-1)^2}];

MatrixForm[h];

{vals,vecs}=Eigensystem[h];

vals[[(n-1)^2]]/π^2
vals[[(n-1)^2-1]]/π^2
vals[[(n-1)^2-2]]/π^2
vals[[(n-1)^2-3]]/π^2

num = (n-1)^2;
vals[[num]]/π^2
data2D = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2D]

NOTE : Doesn' t always pick out eigenmodes with conventional symmetry in case of degenerate states -- depends on choice of 
n (even or odd).   But can combine states to get conventional ones.

num = (n-1)^2-1;
vals[[num]]/π^2
data2Da = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Da]

num = (n-1)^2-2;
vals[[num]]/π^2
data2Db = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Db]

ListPlot3D[data2Da + data2Db]

Plot3D[Sin[2π x]Sin[π y],{x,0,1},{y,0,1}]

NOTE : Eigenvalues good until wavelength approaches grid spacing.

vals/π^2

a = {};
For[i=1,i≤9,++i,For[j=1,j≤9,++j,a = Append[a,i^2+j^2]]]

a;

Sort[a]

vecs[num].vecs[num]

Length[vecs[[1]]]

vecs[[num-4]].vecs[[num-4]]

2-D harmonic oscillator

n = ii = jj = 50;
xl = -4.;
yb = -4.;
xr = 4.;
yt = 4.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
y[i_]:= yb + i Δ
u0 = 4.;

ijToIndex[i_,j_]:= (i-1)*n +j                                                                      (* transform x-y coordinates to 1-D indexing *)
indexToij[i_]:= {Ceiling[i/(n-1)],i-(Ceiling[i/(n-1)]-1)*(n-1)} (* transform 1-D indexing to x-y *)
u[x_,y_]:= u0 (x^2+y^2)

NOTE : Tricky part is getting points on boundary right.

h = Table[If[i==j,4./Δ^2+ u[x[indexToij[i][[1]]],y[indexToij[i][[2]]]],0]
              + If[j==i+1 && Mod[j,n-1]-1 ≠ 0,-1/Δ^2,0]
              + If[j==i-1 && Mod[i,n-1]-1 ≠ 0,-1/Δ^2,0]
              +If[(j== i+(n-1)||j== i-(n-1)),-1/Δ^2,0],{i,1,(n-1)^2},{j,1,(n-1)^2}];

MatrixForm[h];

{vals,vecs}=Eigensystem[h];

vals[[(n-1)^2]]
vals[[(n-1)^2-1]]
vals[[(n-1)^2-2]]
vals[[(n-1)^2-3]]
vals[[(n-1)^2-4]]
vals[[(n-1)^2-5]]

num = (n-1)^2-5;
vals[[num]]
data2D = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2D,PlotRange->All]

NOTE : Doesn' t pick out eigenmodes with conventional symmetry in case of degenerate states.  
But can combine states to get conventional ones.

num = (n-1)^2-1;
vals[[num]]
data2Da = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Da,PlotRange->All]

num = (n-1)^2-2;
vals[[num]]/π^2
data2Db = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Db,PlotRange->All]

ListPlot3D[data2Da + data2Db,PlotRange->All]

2-D harmonic oscillator (change position of Δ’s)

n = ii = jj = 80;
xl = -5.;
yb = -5.;
xr = 5.;
yt = 5.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
y[i_]:= yb + i Δ
u0 = 4.;

ijToIndex[i_,j_]:= (i-1)*n +j                                                                      (* transform x-y coordinates to 1-D indexing *)
indexToij[i_]:= {Ceiling[i/(n-1)],i-(Ceiling[i/(n-1)]-1)*(n-1)} (* transform 1-D indexing to x-y *)
u[x_,y_]:= u0 (x^2+y^2)

NOTE : Tricky part is getting points on boundary right.

h = Table[If[i==j,4.+ Δ^2 u[x[indexToij[i][[1]]],y[indexToij[i][[2]]]],0]
              + If[j==i+1 && Mod[j,n-1]-1 ≠ 0,-1,0]
              + If[j==i-1 && Mod[i,n-1]-1 ≠ 0,-1,0]
              +If[(j== i+(n-1)||j== i-(n-1)),-1,0],{i,1,(n-1)^2},{j,1,(n-1)^2}];

MatrixForm[h];

{vals,vecs}=Eigensystem[h];

vals[[(n-1)^2]]/Δ^2
vals[[(n-1)^2-1]]/Δ^2
vals[[(n-1)^2-2]]/Δ^2
vals[[(n-1)^2-3]]/Δ^2
vals[[(n-1)^2-4]]/Δ^2
vals[[(n-1)^2-5]]/Δ^2

Table[vals[[(n-1)^2-i]]/Δ^2,{i,0,100}]

num = (n-1)^2-5;
vals[[num]]/Δ^2
data2D = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2D,PlotRange->All]

NOTE : Doesn' t pick out eigenmodes with conventional symmetry in case of degenerate states.  
But can combine states to get conventional ones.

num = (n-1)^2-27;
vals[[num]]/Δ^2
data2Da = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Da,PlotRange->All]

num = (n-1)^2-2;
vals[[num]]/Δ^2
data2Db = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2Db,PlotRange->All]

ListPlot3D[data2Da + data2Db,PlotRange->All]

num = (n-1)^2-100;
vals[[num]]/Δ^2
data2D = Table[{x[indexToij[i][[1]]],y[indexToij[i][[2]]],vecs[[num]][[i]]},{i,1,Length[vecs]}];

ListPlot3D[data2D,PlotRange->All]

2 π/Sqrt[12.]

Δ

(2π/Δ)^2

Hydrogen radial equation

Notes : Variable is x, not r, to maintain similar notation to that used in other 1-d problems solved in this notebook. 
            Units for x are a0.
            Solve differential equation for u(r), where ψ(r) = u(r)/r.

xl = 0.;            (* "left" endpoint, r = 0 *)
xr = 4000.;    (* "right" endpoint; must scale with n^2 of highest-energy state desired *)
Δ = 0.5 ;          (* grid-spacing *)
n = xr/Δ           (* number of points *)
x[i_]:= xl + i Δ
l =1.;              (* Orbital angular momentum *)
u[x_]:= -2/x + l(l+1)/x^2   (* Effective potential *)

Fill Matrix:

h = Table[If[i==j,2./Δ^2+u[x[i]],0]+ If[j==i+1 || j==i-1,-1/Δ^2,0],{i,1,n-1},{j,1,n-1}];
MatrixForm[h];

Calculate eigenvalues/eigenvectors:

{vals,vecs}=Eigensystem[h];

ListPlot[vals];

Technique will only work for bound states, so find them:

boundVals = {};
boundVecs = {};
For[i=1,i≤ Length[vals],++i,If[vals[[i]]<0,{boundVals =Append[boundVals ,vals[[i]]],
boundVecs = Append[boundVecs,vecs[[i]]]}]]

boundVals;

Compare calculated values with 1/n^2:

MatrixForm[Table[{-1./(n+l)^2,boundVals[[n]]},{n,l+1,Length[boundVals]}]]

Compare numerical results to analytical results
NOTE: In plots below, I have converted back from u[r] to ψ[r]

a = 1;
rR[n_,l_,r_]:= 
Sqrt[(2/(n a))^3 Factorial[n-l-1]/(2n Factorial[n+l]^3)] Exp[-r/(n a)] (2 r/(n a))^l Factorial[n+l]LaguerreL[n-l-1,2 l + 1,2r/(n a)]
ψ[n_,l_,m_,r_,θ_,ϕ_]:= rR[n,l,r]SphericalHarmonicY[l,m,θ,ϕ]

Results for n = 1, l = 1

nn = Length[boundVals]
num = 1;
data =ListPlot[Table[{x[i],boundVecs[[num]][[i]]/x[i]},{i,1,Length[boundVecs[[num]]]}],PlotRange->{{0,40},All}]

thy =Plot[ rR[num+l,l,r],{r,0,40},PlotRange->All]

Results for n = 41, l = 1, near origin

num = 40;
data =ListPlot[Table[{x[i],boundVecs[[num]][[i]]/x[i]},{i,1,Length[boundVecs[[num]]]}],PlotRange->{{0,40},All}]

thy =Plot[- rR[num+1,1,r],{r,0,40},PlotRange->All]

Results for n = 41, l = 1, near outer turning point

num = 40;
data =ListPlot[Table[{x[i],boundVecs[[num]][[i]]/x[i]},{i,1,Length[boundVecs[[num]]]}],PlotRange->{{1000,xr},{-0.000015,.000015}}]

thy =Plot[- r  rR[num+l,1,r],{r,1000,xr},PlotRange->All]

Helium

See Barraclough and Mooney, J. Chem. Phys., vol. 54, p. 35 (1971); Hajj, J. Phys B, vol. 15, p. 683, (1982)
Differential equation is for ϕ = r1*r2*ψ.

n = ii = jj = 70;
xl = 0;
yb = 0.;
xr = 5.;
yt = 5.;
Δ = (xr-xl)/n;
x[i_]:= xl + i Δ
y[i_]:= yb + i Δ

ijToIndex[i_,j_]:= (i-1)*n +j                                                                      (* transform x-y coordinates to 1-D indexing *)
indexToij[i_]:= {Ceiling[i/(n-1)],i-(Ceiling[i/(n-1)]-1)*(n-1)} (* transform 1-D indexing to x-y *)
u[x_,y_]:= -4/x - 4/y + 2 /Max[x,y]

NOTE : Tricky part is getting points on boundary right.

h = Table[If[i==j,4./Δ^2+ u[x[indexToij[i][[1]]],y[indexToij[i][[2]]]],0]
              + If[j==i+1 && Mod[j,n-1]-1 ≠ 0,-1/Δ^2,0]
              + If[j==i-1 && Mod[i,n-1]-1 ≠ 0,-1/Δ^2,0]
              +If[(j== i+(n-1)||j== i-(n-1)),-1/Δ^2,0],{i,1,(n-1)^2},{j,1,(n-1)^2}];

{vals,vecs}=Eigensystem[h];

ListPlot[vals];

Technique will only work for bound states, so find them.

boundVals = {};
boundVecs = {};
For[i=1,i≤ Length[vals],++i,If[vals[[i]]<0,{boundVals =Append[boundVals ,vals[[i]]],
boundVecs = Append[boundVecs,vecs[[i]]]}]]

boundVals/2

num = 1
boundVals[[num]]/2
data2D = Table[{r1 = x[indexToij[i][[1]]],r2 =y[indexToij[i][[2]]],-boundVecs[[num]][[i]]/(r1 r2)},{i,1,Length[vecs]}];

ListPlot3D[data2D,PlotRange->All]

Convert Hartrees to eV

Look at difference between solution and simple product approximation :

resid= Table[{r1=data2D[[i,1]],r2=data2D[[i,2]],data2D[[i,3]]- data2D[[1,3]]Exp[-1.5(r1-Δ+r2-Δ)]},{i,1,Length[data2D]}];

ListPlot3D[resid,PlotRange->All]