Now I have the output \$v(t)\$ as array of sampled voltages.
It seems like you are missing a fundamental concept about signal processing. An array is not a function of time. A function \$V(t)\$ is dependent on time, but when you have an array of numbers in a software program, time is completely irrelevant. They are just values stored in memory. When you access the values in the array, they become a function of the index. You can call these values samples, and say they are dependent on the sample index n. This means what you really have is \$v[n]\$ where there are n samples, and \$v[n]\$ is the value of the nth sample.
Now how you obtain the samples is a different story. If \$v[n]\$ is a discrete voltage signal, you likely used some method to measure the voltage of a continuous voltage signal \$V(t)\$ every T seconds, and recorded those values. Remember the recorded values are not a function of time, but a function of index. Once the values are recorded, all you have is a list of numbers with no reference to time whatsoever. You can make those values meaningful with respect to time if you remember the sampling interval T. Then you can make the following idealized association: \$v[n] = V(nT)\$. Practically however there is likely some measurement error that makes the above association only almost true.
how can I obtain the input I(t) as array of samples?
Here is one approach:
Convert your transfer function \$H(s)\$ into the equivalent discrete-time system \$H(z)\$. There are many ways to map the s-domain into the z-domain. Some ways are more accurate than others. A pretty standard approach is to use the bilinear transform:
$$s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$$
Find the frequency response of \$H(z)\$ (\$H(z) = H(j\omega_d)\$ where \$\omega_d\$ is the discrete-time frequency of the signal). Typically, you can calculate \$\omega_d\$ as \$\omega_d = \omega T\$, but the bilinear transform causes non-linear frequency warping, so you need to use the following formula instead: $$\omega_d = \frac{2}{T} arctan(\frac{\omega T}{2})$$ When you use this in Python, you will need to represent \$\omega_d\$ with a finite array which makes it a discrete function \$H[k]\$.
Take the DFT of \$v[n]\$ to get \$V[k]\$ which is a discrete representation of the DTFT of \$v[n]\$ (which is \$V(j\omega_d)\$). In Python all you need to do is numpy.fft.fft(v)
where v
is your array of voltages.
Use the deconvolution theorem to calculate \$I[k]\$. You mentioned in the question that you know $$V(j\omega) = I(j\omega) H(j\omega)$$ This means you can simply rearrange so (for the discrete case) $$I[k] = \frac{V[k]}{H[k]}$$ Here's a stack overflow question/answer on how to do deconvolution in Python.
Use the inverse DFT to calculate \$i[n]\$.
Now if you're just looking for the samples that correspond to \$I(t)\$ at every interval \$I(nT)\$, then you can stop there. If you really want to find \$I(t)\$ in continuous time, you will have to reconstruct it from \$i[n]\$.