multiple output array to calculate Normalisation of wave function

Hello, I'm trying to create a function that allows me to check the normalisation of the wave function as is done in Quantum Mechanics (See Griffiths intro to QM)

I start by calculating the constant coefficients, also taking their respecitve reciprocals. These are the constants I will use when checking for normalisation.

c = zeros(1,6);
n = [1 2 3 4 5 6];
for j = 1:length(n)
c(j) = cn(n(j),x);
c(j) = 1./c(j);
end

function out = cn(n,x)
a = max(x);
fun = @(x) sin((n*pi*x)/a).*psi_n(1,x);
 
out = sqrt(2/a)*integral(fun,-a,a);

end

I've just chucked both scripts in to one here, they are separate files though.

Next I use a script to calculate the wave function for some specific n. As this will be squared during the normalisation check, I've left the time dependence out. (in case you're wondering).

function out = psi_n(n,x)
    a = max(x);
    out = sqrt(2/a)*sin((n*pi*x)/a);
end

Now I want to consolidate this into one script and calcluate multiple wave functions at once. So for example I come with an n vector - n = [1 34 54 709]. The script takes it and calculates the constant coefficient, then the whole wave function for each x and outputs them individually. i.e c(n)*ψ(1 ,34, 54, 709). The last part were I integrate and square the wave function I think I can figure out.

I tried the following, and I know it doesn't work because of the assignment values are wrong on each side of the equals sign but I'm not sure how to fix the output...

function out = normalisation_checker(x,c)

n = [1];
a = max(x);
psi = zeros(1,256);
    for j = 1:length(n)
        psi(j) = psi_n(n(j),x);
    out = (c(j)).*sqrt(2/a).*psi_n(n(j),x);
    end
    
end