An algorithm to investigate numerically the solitary wave-like solutionsof Boussinesq differential equation is constructed. The bifurcation problem is re-formulated as an inverse problem for coefficient identification, introducing a newparameter. Such a way the nontrivial solution is separated from the trivial one.The Method of Variational Imbedding (MVI) is used for solving the inverse prob-lem. As illustration of the approach a Mathematica code with the finite differencescheme is prepared and the numerical results, including verification of the scheme,connection between the solution and the coefficients, and capability of Mathematica to paralelize the algorithm are presented.