I know that up sampling by L means padding L-1 zeros between each samples. My question is how it is done practically? Do we use MATLAB or any other software, or there are hardware circuits to do so.
If padding L-1 zeros were the only thing that's done for upsampling, nobody would actually store the end result with zeros that way, because the zeros would be redundant.
In practice, padding with L-1 zeros is always followed by low pass filtering. If that filtering is done with an FIR filter, a dumb implementation would have many coefficients being multiplied by 0. A smart implementation will instead use a polyphase approach where only those coefficients that won't see a guaranteed 0 will be multiplied and summed.
So that partially answers your question: most hardware or software implementations won't even do any padding. They will instead cycle through a set of polyphase filter coefficients, and do something with the result of these successive FIR filter outputs.
Things might be different when doing Matlab or a NumPy implementation: in that case, quick prototyping is often more important than the computational or storage efficiency of the implementation, and rather than using a polyphase filter, people will indeed insert zeros in between samples and convolve the result with the non-polyphase version of the FIR filter.