function [psiEven,psiOdd] = schroInt2(potType,Erel)
% This function tries to find solutions to the Schrodinger equation for
% various potentials.  The type of potential is determined by the
% argument potType.
% case of potType = 1 : is a square well potential
% case of potType = 2 : is a square double well potential
% case of potType = 3 : is a smooth double well potential
% The parameter Erel is the energy relative to the potential well, and 
% should be between 0 and 1.  
% A value of Erel = 0 sets the energy at the bottom of the well,
% and Erel = 1 sets the energy at the top of the well,
% while Erel = 0.5 is half way between the top and bottom.  

figure(10)
psiEven = schroIntEvenOdd(potType,'even',Erel);
title(sprintf('Even trial function; E=%0.12f',Erel))
figure(11)
psiOdd = schroIntEvenOdd(potType,'odd',Erel);
title(sprintf('Odd trial function; E=%0.12f',Erel))
end

function psiFull = schroIntEvenOdd(potType,even_odd,Erel)    
    [U, Emin, Emax, x, dx] = potential(potType);
    N = length(x);
    E = Emin+(Emax-Emin)*Erel;
    kappa2 = U-E;

    psi = NaN(1,N);
    switch even_odd
        case 'even'
            psi(1) = 1;
            psi(2) = psi(1) + 0.5*dx^2*kappa2(1);
        case 'odd'
            psi(1) = 0;
            psi(2) = dx;
    end
            
    for n = 2:(N-1)
        psi(n+1) = kappa2(n)*psi(n)*dx^2 + 2*psi(n) - psi(n-1);
    end
    %pmax = max(abs(psi(x<0.5*a)));
    psiClassical = psi(kappa2<0);
    if ~isempty(psiClassical)
        pmax = max(abs(psiClassical));
        psi = psi/pmax;
    end
        
    %psi = psi*sqrt(2*psi*psi'*dx); % normalize the wave function;
    
    switch even_odd
        case 'even'
            psiFull = [fliplr(psi(2:N)),psi];
        case 'odd'
            psiFull = [-fliplr(psi(2:N)),psi];
    end
    
    xFull = [-fliplr(x(2:N)),x];
    UFull = [fliplr(U(2:N)),U];
    plot(xFull,UFull/max(abs(U)),'.g','markersize',1);
    hold on
    ind = find(UFull>E); % find the classically forbidden region.
    plot(xFull(ind),UFull(ind)/max(abs(U)),'.r','markersize',1)
    plot(xFull,psiFull)
    hold off
    ylim(1.5*[-1,1])
    xlim([min(xFull),max(xFull)])
end

