clear
% the value of a needs to be adjusted to match the condition from the text
% of the problem 
a = 0.001;
% s is the step size for integration
s = 0.5;
delta = 1;
% ener is the initial guess for the energy. it needs to be varied
ener = 0.03;
iter = 0;
% pocetak specifies half of the integration range
pocetak = 100;
% broj is the total number of integration points
broj = pocetak*2/s + 1;
% delta is the convergence criterion; you need to specify it
while abs(delta)>1e-3
    iter = iter + 1;
    br = 0;
    deo = 0.0000000001*ener;
    for E = 9999999999*deo:deo:10000000001*deo;
        br = br + 1;
        q = 0;
% x is the coordinate
        for x = -pocetak:s:pocetak;
            q = q + 1;
% index 1 refers to left-to-right wavefunction and index 2 right-to-left
% wavefunction
            coord1(q) = x;
            coord2(q) = -x;
% V is the potential; it needs to be modified if you want to use this code
% for the part b
            V(q) = a*abs(x);
            G1(q) = (2*V(q) - 2*E);
            G2(q) = G1(q);
% because the problem is symmetric
        end
% here you need to specify boundary conditions for both the wavefunction 1
% and wavefunction 2
        psi1(1) = % asymptotic value of the wavefunction
        psi1(2) = % value of the wvf at the second point. use what you know about asymptotic behavior of the solutions
        y1(1) = (1-((s^2)/12)*G1(1))*psi1(1);
        y1(2) = (1-((s^2)/12)*G1(2))*psi1(2);
        psi2(1) = % similarly as above
        psi2(2) = % similarly as above
        y2(1) = (1-((s^2)/12)*G2(1))*psi2(1);
        y2(2) = (1-((s^2)/12)*G2(2))*psi2(2);
        for n = 2:broj-1;
            y1(n+1) = 2*y1(n) - y1(n-1) + s^2*G1(n)*psi1(n);
            psi1(n+1) = (1-((s^2)/12)*G1(n+1))^(-1)*y1(n+1);
            y2(n+1) = 2*y2(n) - y2(n-1) + s^2*G2(n)*psi2(n);
            psi2(n+1) = (1-((s^2)/12)*G2(n+1))^(-1)*y2(n+1);
        end
        for n = 1:broj;
            psii2(broj + 1 -n) = psi2(n);
        end
 % the next loop searches for the first lobe from the right at the wvf 2
        for n = broj:-1:2;
            if abs(psii2(n))>abs(psii2(n-1))
 % jp1 is the matching point
                jp1 = n;
                fa1 = psi1(jp1);
                fa2 = psii2(jp1);
            break
            end
        end
 % values of both the wvfs are set to 1 at the matching point
        for n = 1:broj;
            psi1(n) = psi1(n)/fa1;
            psi2(n) = psii2(n)/fa2;
            y1(n) = (1-((s^2)/12)*G1(n))*psi1(n);
            y2(n) = (1-((s^2)/12)*G2(n))*psi2(n);
        end
 % ff is the correction function
        ff(br) = (-y2(jp1+1) + 2*y1(jp1) - y1(jp1-1))*s^(-2) + G1(jp1)*psi1(jp1);
    end
 % izv is the derivative of ff wrt to energy
    izv = (ff(3) - ff(1))/(2*deo);
    delta = -ff(2)/izv;
 % energy gets corrected here
    ener = ener + delta;
end
nin = 0;
plp = ((pocetak - (1.8*ener/a))/s) - mod(((pocetak - (1.8*ener/a))/s),1);
 for ninx = plp:(broj-plp);
     nin = nin + 1;
     psinov1(nin) = psi1(ninx);
     psinov2(nin) = psi2(ninx);
     coordn1(nin) = coord1(ninx);
 end
 % normalization
 pss1 = psinov1.^2;
 pss2 = psinov2.^2;
 psisum1 = sqrt(s*sum(pss1));
 psinov1 = psinov1./psisum1;
 psisum2 = sqrt(s*sum(pss2));
 psinov2 = psinov2./psisum2; 
 plot(coordn1,psinov2)
 hold on
 plot(coordn1,psinov1)
 format long
 ener