How do I assign a double to a symbolic variable given that the scope of the variable does not span the function in which I work?

I have created a function which calculates the Hermite polynomials given two initial values: H0=1 and H1=2x, where x is a symbolic variable. I define a symbolic array initially filled with random symbolic expressions inside this function. Then I update the cells of it using the recurrence relation for Hermite polynomials. Note that H0 corresponds to first cell, H1 to second cell and so on. This function returns a vector of symbolic expressions which I want to use in the main function. I invoke (inside the main function) the function we talked about above and I get that vector of expressions. Then I want to focus on a particular cell of it. I want to update that cell (which contains a symbolic expression) with different values of the symbolic variable. (to simulate the real axis). I use a for loop to do this. However, after every iteration, the result appearing when I call disp(hermitevector(n+1)) is the same, irrespective to the iteration process the index has gone through. I have defined that symbolic expression in the function which generates Hermite polynomials. Does the main function understand about it? Does it know it? If not, I thought I may define it as a global variable using global x. Nothing has changed.

Please see my attached code and give me some indication on how to manipulate symbolic expressions. Thank you very much!!!

 % function [ hermitevector ] = hermite( uppervalue )  %want x to be symbolic. put an uppervalue to have a broad range of Hermite polynomials, up until the required n. also attach n such that we can link this input para to the actual n in the main script
% hermite function iteratively calculates the values of the Hermite
% polynomials starting from two intial values: H0=1 and H1=2x.
x = sym('x');
hermitevector = sym('hermitevector',[1,uppervalue+1]);  %creates a vector smth like: [hermitevector1, hermitevector2, hermitevector3,... hermitevectoruppervalue];
%global x;
    hermitevector(2)=2*x;
    hermitevector(1)=1;
for j=1:(uppervalue-1)  %create a vector of symbolic expressions, as x is declared to be symbolic. it runs up until it gets the Hermite_uppervalue (this should be n).
    hermitevector(j+2)=2*x*hermitevector(j+1) -2*j*hermitevector(j);
end
end

main function I was talking about:

function [ psitheo, psi ] = main_script(delta, x1, n, E)
%this function integrates numerical Schrodinger equation from lower value
%x0 to upper value x1, plotting both the theoretical result and the
%computed result.
f = @(x) (x.*x-E);   %define what is f. basically the RHS of the numerical Schrodinger equation. potential=x^2 . energy=E 
%&&&&analytical solution&&&&
if ( mod(n,2)==0 )
    psi0=1; dpsi0=0;  %set boundary conditions
else
    psi0=0; dpsi0=1;  %set boundary conditions
end
 psitheo = sym('psitheo',[1,n+1]); %declare the psitheoretical as a row vector of dimension x1, initially filled random symbolic expressions
 psi= sym('psi',[1,n+1]); %declare the psi as a row vector of dimension x1, initially filled with random symbolic expressions

%arrayofintegration=0:0.01:x1; %simulate real axis 
hermitevector = sym('hermitevector',[1,n]); %create a symbolic vector full of symbolic expressions for the Hermite polynomials
hermitevector = hermite(n+1);  %populate it with Hermite polynomials of x (x symbolic variable). input para:n --> need n-th order Hermite poly, so acces n+1 cell of hermitevector[] array, as indexing starts at 1 and first cell holds H0.
noofsteps=x1/0.01; %real axis simulated here
track= 0:0.01:x1;  %for later purposes: graph plotting
for k=1:(noofsteps+1) %densely packed values of x, ranging from 0 to x1
    hermitevector(n+1)=subs(hermitevector(n+1),(k-1)*delta);%PROBLEM! always and always hermitevector(n) is a specific value, for any k. it should change.update the n-th Hermite poly by replacing the symbolic variable x with a real number denoted--> this should yield a double. it is the same as in the script where x=j*delta. their j is our k
    disp(hermitevector(n+1)); %should give different values as k increases.
    psitheo(k)=hermitevector(n+1)*exp(-(k*0.01)^2/2); %fill the vector of psitheoretical. indexing uses k. however, x=k*delta. use the form of solution given in script. hermite runs 'in background' and computes the nth Hermite poly. n given as input variable in this script
end
%disp(size(track)); disp(size(psitheo));
figure
    subplot(2,1,1)
    plot(track,psitheo); 
    title('Theoretical Psi for given n'); 
    xlabel('x'); 
    ylabel('Psi');

%&&&&numerical solution&&&&
%if E==2*n+1   %if we have a closed form solution
 %   psi = solve_numerov(f,track,psi0,dpsi0);
%end

%disp('size of psi vector:'); 
%disp(size(psi)); sizepsi=size(psi);
%sizeoftrack=size(track); 
%disp('track vector'); disp(track);
%final=sizepsi(1,2);
%disp('final'); disp(final);

%for i=1:final
   %psi(i)=subs(psi(i),delta);%update the value of delta we have as input for every cell in the psi[] array. however, this gives all psi(i) equal to 0.05
   %disp(psi(i));
%end
%disp('sizeoftrack(1,2)');disp(sizeoftrack(1,2));
%disp('type of psi'); disp(class(psi(200)));
%subplot(2,1,2);
%plot(track,psi); title('Numerical simulation for Psi'); xlabel('x'); ylabel('Psi');

% end
end

If anyone has the patience to go through the code(which looks big, but it is not complicated) I would much appreciate to tell me what I am doing wrong in the theoretical part, as the function solve_numerov.m is:

function [ psi ] = solve_numerov( f, x, psi0, dpsi0 )
% Solving d ˆ 2 ( psi )/dx ˆ2 = f ( x ) psi from x0 to x1 with boundary
% conditions:
% psi ( x0 ) = psi0 and d( psi ( x0 ) )/dx = dpsi0;
% Input:
% * f: The function to be called on the right hand side of the equation
% (see hint section ), receives x and returns the value f(x).
% * x: array of integration, the first element is the lower bound, x0, and
% the last element is the upper bound, x1.
% * psi0: the value of \psi at x0, i.e. \psi(x0)
% * dpsi0: the value of \psi derivative at x0, i.e. d\psi(x0)/dx
%
%Output:  
% * psi: array of values of \psi with each element in the array
% corresponding to the same element in x

%Example use:
%>> f = @(x) x^2-1;
%>> x = linspace(0, 1, 100);
%>> psi0 = 1;
%>> dpsi0 = 0;
%>> psi=solve_numerov(f, x, psi0, dpsi0);
%>> plot(x, psi);
%f = @(x) x.*x-E;%define an annonymous function to calculate the value of the RHS term in Schrodinger equation. E is constant and is input para for main_script.m
sizex=size(x);%sample the array of integration such that we are able to write the argument for psi function at every step of iteration
psi=sym('psi',[1,sizex(1,2)]);  %declare this vector array as a symbolic one.
deltaa = sym('deltaa');  %do not know what is deltaa now. deltaa is input parameter in main_script.m, so leave it as symbolic variable and update it in main_script.m
%global deltaa;
%psi(1)=psi0; %put it in the vector array
if (psi0==1)  %check parity of solutions
psi(2)= psi0 + (deltaa^2/2)*f(0)*psi0 + (deltaa^4/24)*(f(0)^2*psi0); %for even n(n is input parameter in main script)
elseif (psi0==1)
psi(2)= deltaa*dpsi0 + (deltaa^3/3)*f(0)*dpsi0; %for odd n(n is imput parameter in main script)
end

for k0= 1:(x-1)
    psi((k0+1))=((2 +(5/6)*deltaa^2*f(k0*deltaa))* psi(k0) - (1-((deltaa^2)/12)*f((k0-1)*(deltaa)))*psi(k0-1)) /(1 -(deltaa^2/12)*f((k0+1)*(deltaa)));
end    

end