creating a bcd squarer using verilog
Basically Verilog operations are done in binary not BCD. So if you do something like i**2 the result will be in binary. Splitting the nibbles doesn't convert binary to BCD.

You have two options. Either:


Continue and implement a binary to BCD conversion - you could do this as a function, which because it is used in an initial block shouldn't infer any hardware assuming it is something the synthesizer can calculate. Or
Compute all the values before hand, and simply initialise the LUT you have designed with the precomputed values. This is by far the easiest way, and to be honest likely (or so I've found) the usual way of doing it. 


In the latter option you can do this several ways. One is to use something like Excel to compute a table of values and initialise the LUT from that - e.g. by reading the values in turn from a file (Verilog has macros for this like $readmem).

Alternatively if you are feeling adventurous, you can create a script (e.g. TCL) which creates the values and generates either the a Verilog file containing the code for the whole LUT including initialation, or creates a memory initialisation file which you can load into your LUT. This option is not really useful for this particular example, but great for more complex ones, perhaps save it for once you've gained more experience.