How to use numerical methods to determine the eigenvalues of TISE?

I've successfully used numerical methods to derive the TISE solution for a Simple Harmonic Oscillator potential problem. I'm able to plot graphs of the Wave equation.

Is there a way which I can make use of numerical methods to find the eigenvalues of TISE?

Sure, it's already contained in your data. If you already have $\psi_n(x)$, then the relation $$\frac{\hat H\psi_n(x)}{\psi_n(x)} = \frac{-\frac{\hbar^2}{2m}\psi_n''(x) + \frac12m\omega^2 x^2\psi_n(x)}{\psi_n(x)} = E_n$$ gives you the eigenvalue in terms of quantities that you can calculate explicitly. (For better accuracy, I'd discard the lowest quartile in $|\psi_n(x)|^2$, or you're dividing by a small number which can be inaccurate, and then average over the rest.)

Now, depending on what method you actually used, there's probably a better (more accurate, more robust, etc) way to get the eigenvalues as part of the process, but the details of this will depend on the method and it's hard to say anything else that will hold in general.

If you have the wavefunctions you automatically have the eigenvalues since $$-(1/2) \psi''[x] + (1/2) x^2 \psi[x] = \epsilon \psi[x] \tag{1}$$ will converge only when $\epsilon$ is an eigenenergy.

If you do not have $\epsilon_n$, there are several methods available, but all are ultimately based on the observation that boundary conditions are critical in producing quantized eigenvalues of energy.

In the case of the harmonic oscillator, the simplest methods uses the boundary condition is $\psi(\pm \infty)=0$. In principle, only if the energy is exactly right can one satisfy this condition, for otherwise the solution will diverge exponentially at large $x$.

When solving analytically the Schrodinger equation, the sequence is to first the search for solutions that satisfy this condition and then condition to find condition on the energies.

In the simplest numerical schemes the sequence is reversed: one "guesses" an energy, solves the resulting differential equation, and checks to see if the resulting solutions satisfies the boundary solution for sufficiently large $\vert x\vert$. There is always some roundoff error, and so the solution will never be strictly $0$ for very large $x$, but practically it is enough to require the solution to go to $0$ over a reasonable range of $x$ for sufficiently large $x$ before the solution blows up again.

Exactly what "very large" means, and exactly what is a "reasonable range" is a matter of numerical accuracy and may depend on the specifics of the integration schemes. A Runge-Kutta integrator of order $8$ will typically do better than one of order $4$, in the sense that the solution from RK8 will be remain small for longer at large $x$ than the solution from RK4.

In the case of the harmonic oscillator (or any other symmetric potential), the initial conditions on the derivative and the value at $x=0$ are easy: for the even solutions one can take $\psi(0)=1$ and $\psi'(0)=0$ while for the odd one $\psi(0)=0$ and $\psi'(0)=1$ are appropriate. These initial conditions do not affect the accuracy of the solution. The solution is not normalized but the normalization is unimportant to accurately determine the energy.

Thus given an initial guess $e_0$, one launches the integrator and checks to see if the solutions blows up for values of $x$ "not large enough". If the solution blows up one makes another guess, say $\epsilon_1> \epsilon_0$. If the new solution is worse, then redo with $\epsilon_1<\epsilon_0$, and this way by trial and error until a guess value is found so that the solution remains small over the desired range.

This is illustrated in the two figures below, which are numerical solutions to $$-(1/2) \psi''[x] + (1/2) x^2 \psi[x] = (n + 1/2) \psi[x]$$ with initial conditions $\psi=1$ and $\psi'=0$. The graph on the left used a Runge-Kutta integrator of order 8, and the one on the right a RK scheme of order 4. In both graph, the black plots are for the exact value $\epsilon=2$ (in units of $\hbar\omega$), the red plots for $\epsilon=1.95$ and the blue plots for $\epsilon=2.05$. This clearly shows how the integration is quite sensitive to the initial guess energy, how the solution "rapidly" diverges when the guess energy is wrong, and how the asymptotic tail of the solution improves with the order of the integrating scheme. (You can see a hint, on the left, that even the exact solution $\epsilon=2$ produces a solutions that is about to diverge beyond $\vert x\vert=7$.)

Because the number of nodes (or zeroes) of each eigenfunction increases with energy, it is easy to verify that no energy value has been missed. The lowest energy solution will have no node, the first excited state one node etc. Thus, if a guess energy produces a valid solution with more nodes than expected, one simply tries a smaller guess. In this way, and with patience, one can find with reasonable accuracy all energy eigenvalues, especially if an accurate integrator is used.

There are a variety of other more sophisticated schemes. One is called the shooting method. Another common method is called the matching method: instead of starting at the centre and checking the solution at large $x$, starts with large $x$ and works backwards towards the centre. When done starting at $+x$ and $-x$, the continuity of the solution and of its derivative is used to evaluate if the guess energy is to be kept or rejected.