Solve the finite potential well problem using Mathematica.
 The solution is based on the following text book:  Mathematica for Theoretical Physics, Gerd Baumann.

(* Solve the Schroedinger equation in regions 1 and 3 x<-a and x>a. *) 
(* The solution is stored in s13. The constant κ is defined as κ = (Sqrt[2m(V0-E)])/ℏ *)
s13=DSolve[(∂_(x,x))ψ[x]-(κ^2)ψ[x]==0,ψ,x]//Flatten


(* Solve the Schroedinger equation in region 2: -a<x<a. *)  
(* The solution is stored in s2. The constant k is defined as k=(Sqrt[2mE])/ℏ *)
s2=DSolve[(∂_(x,x))ψ[x]+(k^2)ψ[x]==0,ψ,x]//Flatten

(* Construct the wavefunctions and define the appropriate integration constants in the various regions. *)

(* In region 1, we must ensure that the wavefunction remains finite for all x and as a result C[2] must be zero. *) 
ps1=ψ[x]/.s13/.{C[1]->A1,C[2]->B1}/.B1->0


(* In region 3, we must ensure that the wavefunction remains finite for all x and as a result C[1] must be zero.  *) 
ps3=ψ[x]/.s13/.{C[1]->A3,C[2]->B3}/.A3->0


(* In region 2, there are no constraints on the integration constants. *)
ps2=ψ[x]/.s2/.{C[1]->A2,C[2]->B2}

(* A valid solution to the Schroedinger equation requires that the wavefunctions and their first derivatives match at x = -a and at x = a. *)

(* This matching condition will generate 4 equations that must be satisfied: eq1, eq2, eq3, and eq4. *)

(* Equation 1: match wavefunctions at x = -a. *)
eq1=ps1==ps2/.x->-a

(* Equation 2: match derivative of wavefunctions at x = -a. *)
eq2=(∂_x)ps1==(∂_x)ps2/.x->-a

(* Equation 3: match wavefunctions at x = a. *)
eq3=ps3==ps2/.x->a

(* Equation 4: match derivative of wavefunctions at x = a. *)
eq4=(∂_x)ps3==(∂_x)ps2/.x->a

(* These four equations can be combined and written in the following way:

(({{., ., ., .}, {., ., ., .}, {., ., ., .}, {., ., ., .}}))(({{A1}, {A2}, {B2}, {B3}}))=0 

There is a non-trivial solution to this equation if the determinant of the matrix is equal to 0.
*)

(* Construct the determinant of the matrix. *)
det1=Map[Coefficient[({eq1,eq2,eq3,eq4}/.Equal[a_,b_]:>a-b),#]&,{A1,A2,B2,B3}]//Transpose//Det//Simplify

(* Determine when the determinant is equal to 0.
The solutions are stored in the spectral array. *)
spectral=MapAll[PowerExpand[#]&,Simplify[Flatten[Solve[det1==0,κ]]]]//FullSimplify

(* Determine the coefficients of the wavefunction for the first solution. *)
sol1=Solve[{eq1,eq2,eq3,eq4}/.spectral[[1]],{A1,B2,A2,B3}]//Simplify//Flatten

(* Determine the coefficients of the wavefunction for the second solution. *)
sol2=Solve[{eq1,eq2,eq3,eq4}/.spectral[[2]],{A1,A2,B2,B3}]//Simplify//Flatten

(* In each region we obtain two solutions: 
one symmetric and one asyymetric solution. 
*)

(* Start with region 1: x < -a *)
ψ1a=ps1/.sol1/.spectral[[1]]
ψ1s=ps1/.sol2/.spectral[[2]]

(* Continue with region 2: -a < x < a *)
ψ2a=ps2/.sol1/.spectral[[1]]
ψ2s=ps2/.sol2/.spectral[[2]]

(* Finish with region 3: a < x *)
ψ3a=ps3/.sol1/.spectral[[1]]
ψ3s=ps3/.sol2/.spectral[[2]]