The reality here is that the FT and DFT are actually two different things. They each have their own consistent theory, but the DFT can only ever be an approximation of the FT. Because of sampling and aliasing, it is easy to find simple cases where the DFT is not "correct" when compared to the FT.
Taking a step back from the math, the mapping between your IFT and IDFT is:
Continuous time $t$ in seconds $\leftrightarrow$ Discrete time sample index $m$ in samples
Continuous radian frequency $\omega$ in rad/sec $\leftrightarrow$ Discrete radian frequency bin center at $\dfrac{2\pi k}{n}$ in rad/sample
To convert your discrete frequency bin index, $k$, into a continuous time frequecy, $f = \dfrac{\omega}{2\pi}$, you'll use the sampling frequency, $f_s = \dfrac{1}{\Delta t}$, thusly:
$$\dfrac{\omega}{2\pi} = f = \dfrac{\frac{f_s}{2}}{\pi}\cdot \dfrac{2\pi k}{n}= \dfrac{f_s}{n}k = \dfrac{k}{n\Delta t}\quad \mathrm{for}\space \dfrac{-n}{2} \le k < \dfrac{n}{2}$$
The range for $k$ highlights one way in which the DFT is not the FT: the DFT is cyclic. The frequency indices $k \ge \dfrac{n}{2}$ represent negative frequencies.
Computing the IDFT also results in cyclic time domain data. If you were planning on computing consecutive IDFTs and concatenating them back together in the time domain, you have to be very careful to get correct results. E.g. For "fast convolution", using FFTs to perform filtering, the Overlap-Add and Overlap-Save algorithms were developed, to deal with aliasing effects, when processing a continual stream of samples when using the block-oriented FFT & IFFT.
Update in response to the questions in the comment
Before I address the questions I will note, that in your original question, you used two different index variables for your discrete time domain sample index: $m$, by your use of $a_m$, and $j$, by your use of $t_j$. I'm going to drop the $j$ and just use $m$.
"... why sampling frequency is $f_s = \dfrac{1}{\Delta t}$ ..."
It is from the definition of frequency and the period at which your time domain is sampled.
In your approach section, for discrete time, you used
$$t_m = t_0 + m\Delta t$$
which, when applied to an input function $f(t)$, yields:
$$f(t_m) = f(t_0 + m\Delta t) = f[m]$$
So $\Delta t$ is the period between the periodic sampling events that yield the sample sequence $f[m]$. $\Delta t$ is your sample period. The frequency of any periodic event is defined as the reciprocal of the period of that periodic event. Therefore
$$f_s = \dfrac{1}{\Delta t}$$
I'll also note for later, that from Nyquist's results, for sampled signals, the frequency range that can be discretized unambiguously is $[f_{min}, f_{max}) = \left[-\dfrac{f_s}{2}, \dfrac{f_s}{2}\right)$, before aliasing back on to frequencies in that interval.
"... why $\omega_k$ should be discretized as $\omega_k = \dfrac{2\pi k}{n} ...$"
I'll present two ways to view why this decision is meaningful.
The first way is to simply compare the kernels of the two different transforms and make a 1-to-1 mapping of the components in the kernels.
A(n Inverse) Fourier Transform kernel is composed of an exponential with the following argument structure:
$$ \pm i (\mbox{time domain variable}) (\mbox{frequency domain variable}) (\mbox{unit conversion constants})$$
When making a comparison to map a IFT kernel to and IDFT kernel:
$$e^{i\omega t} \leftrightarrow e^{i \frac{2\pi k m}{n}}$$
we already have a clearly stated relationship between continuous time $t$ and the discrete time sample index ${m}$, so it is simple enough to keep that as a 1-to-1 mapping.
The remaining continuous radian frequency variable, $\omega$, would then map to the remaining discrete time variable, $k$, and the units conversion constant $\dfrac{2\pi}{n}$ with an additional $f_s$ unit conversion constant in there to undo the normalization.
That tells you how to make a mapping between two different transforms working on two different types of input. However, I suppose you might not find it satisfying, as it doesn't explain why that may be a preferred mapping.
The second way I'll show why one might want to discretize $\omega$ as $\omega_k = \dfrac{2\pi k}{n}$, is to state (without proof) how the transforms manifest in their transform domain complex planes, and then try to show why the $\omega \leftrightarrow \dfrac{2\pi k}{n}$ mapping makes sense.
The Laplace Transform can be thought of as a continuation of the Fourier Transform to (at least the RHS) of the complex s-plane. The Fourier Transform itself is located on the imaginary axis of the s-plane, with $\omega$ normally used as the variable for the ordinate (imaginary) axis.
The DFT can be thought of as being located on the unit circle of the finite sequence Z Transform's z-plane; the unit circle being the curve $z = e^{i\Omega}$. We'll call $\Omega$ the "normalized digital radian frequency", which can only unambiguously range over an interval of length $2 \pi$, e.g. $[0, 2\pi)$ or $[-\pi, \pi)$, before wrapping around and aliasing onto other values of $\Omega$.
So knowing that, for a sampled signal, we can only unambiguously determine frequencies in the range $\left[-\dfrac{f_s}{2}, \dfrac{f_2}{2}\right)$, and we can only have an unambiguous $\Omega$ in a range of $2\pi$, e.g. $[-\pi, \pi)$, we can use the following mapping from s-plane to z-plane to map $\omega$ in $\left\{-2\pi \dfrac{f_s}{2}, 2\pi\dfrac{f_s}{2}\right\}$ to $\Omega$ in $[-\pi, \pi)$:
$$z = e^{s\Delta t} = e^{\frac{s}{f_s}}$$
and for $s = i\omega$, where the Fourier Transform resides in the s-plane, this mapping becomes:
$$z = e^{i\omega\Delta t} = e^{i\frac{\omega}{f_s}}= e^{i\frac{2\pi f}{f_s}}$$
We now have a mapping from the s-plane to z-plane for continuous frequency values that we can unambiguously resolve from a sampled signal. Take note of the ratio $\dfrac{f}{f_s}$ and remember that the unambiguous range of $f$ can only be of length $f_s$, e.g. $\left[-\dfrac{f_s}{2}, \dfrac{f_s}{2}\right)$, before wrapping around and aliasing onto other values of $f$.
So now let's look at the values of the quantity $\dfrac{2\pi}{n}k$. I hope you can recognize them as the $n^{th}$ roots of unity. Those $n^{th}$ roots of unity can range unambiguously for an index set $k$ of length $n$, e.g. $k \in \{0, \dots, n-1\}$ or $k \in \left\{-\dfrac{n}{2}, \dots, \dfrac{n}{2}-1\right\}$, as using a larger set results in wrapping around and aliasing onto other $n^{th}$ roots of unity.
So the quantity $e^{i2\pi\frac{f}{f_s}}$ represents the normalized continuous radian frequencies we can resolve unambiguously when mapped to the z-plane's unit circle continuously. If one has to map that to the discrete values of the IDFT's kernel $e^{i2\pi\frac{mk}{n}}$, on the z-plane's unit circle, one could use the knowledge that the ratios $\dfrac{f}{f_s}$ and $\dfrac{k}{n}$ must both only range over an interval of length 1 to be unambiguous and to avoid aliasing. This naturally leads one to map:
$$\dfrac{2\pi f}{f_s} \leftrightarrow \dfrac{2\pi k}{n}$$
Best Answer
Suppose $u(t, x)$ a solution to the original pde, then define $v(t, x)= e^{-ct}u(t, x)$. Observe \begin{align} v_t-k v_{xx} =&\ -ce^{-ct}u(t, x)+e^{ct}u_t(t, x)-ke^{ct}u_{xx}(t, x) \\ =&\ e^{ct}(u_t(t, x)-ku_{xx}(t, x)-cu(t, x)) = 0 \end{align} which is just the heat equation. Now solve it with Fourier transform or whatever you like. Then we see that \begin{align} u(t, x) = e^{ct}v(t, x) \end{align} is your solution.
Edit: Consider \begin{align} w(t, x) = \exp\left(\frac{c^2t}{4k}+\frac{cx}{2k}\right) u(t, x) \end{align} then we see that \begin{align} \partial_tw-k\partial_{xx}w =&\, \frac{c^2}{4k}\exp\left(\frac{c^2t}{4k}+\frac{cx}{2k}\right) u+\exp\left(\frac{c^2t}{4k}+\frac{cx}{2k}\right)\partial_t u - \exp\left(\frac{c^2t}{4k}+\frac{cx}{2k}\right)\left(\frac{c^2}{4k} u+cu_x+ku_{xx}\right)\\ =&\, \exp\left(\frac{c^2t}{4k}+\frac{cx}{2k}\right)\left(\partial_t u-c u_x-ku_{xx} \right) = 0. \end{align}