Reshaping of 3-dimensional array to construct an isosurface

I'm having an issue finalizing the code that I've ran past line 40 (code below). After reshaping my matrices to perform operations, I am unable to reshape back into an appropriate size to create an isosurface. I'm thinking I may need to slice, or perhaps I'm glossing over something entirely simple. I'm also having issues running the eig function under GPU execution for computational efficiency, but that's tertiary to the problem at hand. If you wish to execute the code to get an idea, I would suggest lowering the N value, as I have, to below 50.

Code:

%establishing iteration size
N = 25;
%establishing meshgrid spacing of appropriate bohr radius to avoid infinite
%square well
x = linspace(-25,25,N);
y = linspace(-25,25,N);
z = linspace(-25,25,N);
[X,Y,Z] = meshgrid(x,y,z);
dx = diff(x(1:2));
%estimating potential of hydrogen atom
e = 1*10^(-10);
d = (sqrt(X.^2 + Y.^2 + Z.^2 + e)).^-1;
n = -dx^2;  
V = d.*n;
%second partial using kronecker product
D = ones(N,1); 
D = [D, -2*D, D];
A = spdiags(D,-1:1,N,N);
%A = full(A);
SPA = sparse(kron(A,A));
%kinetic operator
T = -.5*sparse(kron(SPA,A));
V = reshape(V,1,N^3);
U = sparse(diag(V));
%hamiltonian
H = T + U;
%%H = gpuArray(H);
%calculating discrete eigenstates(energies) in three dimensions, assuming
%no progression in time
[L,C,R] = eigs(H,25,'smallestreal');
%where I'm running into issues
L = reshape(L,1,N^3);
isosurface(x,y,z,L);

Answer 1

[Bruno Luong]

I just fix the code without knowing what is really your aim. I'm surprise you do all kind of kron and eigs and get lost in the underlined dimensions of the result.

I also remove the unnecessary conversion back and forth between full ans sparse.

%establishing iteration size
N = 25;
%establishing meshgrid spacing of appropriate bohr radius to avoid infinite
%square well
x = linspace(-25,25,N);
y = linspace(-25,25,N);
z = linspace(-25,25,N);
[X,Y,Z] = meshgrid(x,y,z);
dx = diff(x(1:2));
%estimating potential of hydrogen atom
e = 1*10^(-10);
d = (sqrt(X.^2 + Y.^2 + Z.^2 + e)).^-1;
n = -dx^2;  
V = d.*n;
%second partial using kronecker product
D = ones(N,1); 
D = [D, -2*D, D];
A = spdiags(D,-1:1,N,N);
%A = full(A);
SPA = kron(A,A);
%kinetic operator
T = -.5*kron(SPA,A);
V = reshape(V,1,N^3);
U = diag(V);
%hamiltonian
H = T + U;
%%H = gpuArray(H);
%calculating discrete eigenstates(energies) in three dimensions, assuming
%no progression in time
NbEig = 4; % 25 for you
[L,C,R] = eigs(H,NbEig,'smallestreal');
for k=1:NbEig
    subplot(2,2,k)
    Lk = reshape(L(:,k),N,N,N);
    isosurface(x,y,z,Lk);
end

[Bruno Luong]

Supposingly the expected figure is similar to this

[Ethan Wilson]

Yes, I believe your method retrieves the appropriate cross section, but there is an error within the operators. Let me look further. I've also realized your axis scaling varies, which is interesting.