Why can't I plot my two-variable symbolic expression using fsurf?

I'm attempting to solve and plot a 1D Time-Dependent Schrodinger equation (TDSE) for an infinite square well symbolically using MATLAB. (I've been away from quantum mechanics and MATLAB for years, so please forgive my ignorance.) I am able to solve the equation correctly, but I can't seem to get fsurf to work. Here is my code:

 clc; close all;
syms x m a hbar pi n A t
n = 2;
psi = sym(zeros(n,1)); 
c = sym(zeros(n,n));
E = sym(zeros(n,1));
Psi = sym(zeros(n,1));

 for i = 1:n
    psi(i) = sqrt(2/a)*sin(i*pi*x/a); %infinite square well wavefunction
end

 Psi0 = A*psi;
A = abs(solve(int(Psi0'*Psi0,x,[0 a])==1,A)); %solves for normalization constant
A = double(A(1)); %changes the symbolic A into a numerical A
Psi0 = A*psi; %redefines the initial wavefunction w/numeric A

 %The below loop acts as the kronecker delta function
for i = 1:n
    E(i) = (n^2*pi^2*hbar^2)/(2*m*a^2); %possible energy states
    for j = 1:n
        if i~=j
           c(i,j) = 0; 
        else
           c(i,j) = psi(i)*Psi0(j);
        end
    end
end

 c = sym(int(sum(c),x,[0 a]))'; %spits out the expansion coefficients

 for i = 1:n
    Psi(i) = c(i)*psi(i)*exp(-1i*E(i)/hbar*t); 
end

 Psi = sym(sum(Psi)); %the most general solution to the TDSE

 %values for substitution
m_elec = 9.10938356*10^-31; %mass of electron
a_ = 10^-10; %angstrom
hbar_ = 1.0545718*10^-34;
pi_ = 3.14159265359;
t_ = [0:0.1:100]';
x_ = [0:a_/(10^3):a_]';
Psi = subs(Psi,{m a hbar pi},{m_elec a_ hbar_ pi_});
% Psi = vectorize(Psi)

 figure('color','white')
fsurf(@(x,t) Psi,[0 2*a_ 0 2*0 a_])

The error I get is as follows:

 Warning: Function behaves unexpectedly on array
inputs. To improve performance, properly vectorize
your function to return an output with the same
size and shape as the input arguments. 
> In matlab.graphics.function.FunctionSurface>getFunction
  In matlab.graphics.function.FunctionSurface/updateFunction
  In matlab.graphics.function.FunctionSurface/set.Function
  In matlab.graphics.function.FunctionSurface
  In fsurf>singleFsurf (line 261)
  In fsurf>@(f)singleFsurf(cax,{f},extraOpts,args) (line 227)
  In fsurf>vectorizeFsurf (line 227)
  In fsurf (line 200)
  In Infinite_Square_Well (line 53) 
 Warning: Error updating FunctionSurface.

 The following error was reported evaluating the
 function in FunctionLine update: Unable to
 convert expression into double array.

I've been working on this for hours, but I can't seem to figure out what to do. I want to generate a 3D plot of Psi(x,t). Also, if you know of a way to keep all my symbols (such as a, hbar, pi, etc) in the equation without substituting numbers in, I would appreciate help with that too.

Answer 1

[Walter Roberson]

f = matlabFunction(Psi, 'vars', [x t]);
fsurf(f, [0 2*a_ 2*0 a_])

However, this will produce an error about index exceeding matrix dimensions. I suspect that has to do with the quite small range of values you are plotting only, in the range of 10^-10 . If you increase that range then the error about matrix dimensions goes away.

The plot that is produced will be empty. That is because the values of the function are complex valued and fsurf() only plots real-valued locations

If you

f = matlabFunction(imag(Psi), 'vars', [x t]);
fsurf(f, [0 1e-8 0 1e-8])

then unfortunately that takes rather long. You would probably be better off using meshgrid() or ndgrid() to create a mesh of locations and evaluate f at those locations.

[Muoio117]

Thanks for your help. I'm not quite familiar with meshgrid(), especially for symbolic functions, but I will look into finding a way to plot my function with it.