There are already several reasonable answers, but to try and add some clarity, here I will try to synthesize these, and do some verification on the resulting algorithm.
Homogeneous Process
For a homogeneous 1D Poisson process with constant rate $\lambda$, the number of events $n$ in a time interval $T$ will follow a Poisson distribution
$$n\sim\text{Poiss}_\lambda$$
And the arrival times $t_1,\ldots,t_n$ will be i.i.d. uniformly through the interval
$$t_i\sim\text{Unif}_{[0,T]}$$
while the waiting times between events will be i.i.d. exponentially
$$\tau_i=t_{i+1}-t_i\sim\text{Exp}_\lambda$$
with CDF
$$\Pr\big[\tau<\Delta{t}\big]=1-\exp^{-\lambda\Delta{t}}$$
So a homogeneous Poisson process can be easily simulated by first sampling $n$ and then sampling $t_{1:n}$ (or, alternatively sampling $\tau$ until $t=\sum\tau>T$).
Inhomogeneous Process
For an inhomogeneous Poisson process with rate parameter $\lambda(t)$ the above can be generalized by working in the transformed domain
$$\Lambda(t)=\int_0^t\lambda(s)ds$$
where $\Lambda(t)$ is the expected cumulative number of events up to time $t$. (As noted in Aksakal's answer, and the reference cited in jth's answer.)
In the generalized approach we simulate in the $\Lambda$ space, where "time" is dilated so in the deformed timeline the Poisson process is now a homogeneous process with unit rate. After simulating this $\Lambda$ process, we map the samples $\Lambda_{1:n}$ into arrival times $t_{1:n}$ by inverting the (monotonic) $\Lambda(t)$ mapping.
Note that for the piecewise-constant rates $\lambda(t)$ here, the mapping $\Lambda(t)$ is piecewise linear, so very easy to invert.
Example Simulation
I made a short code in Octave to demonstrate this (listed at the end of this answer). To try and clear up questions about the validity of the approach for small rates, I simulate an ensemble of concatenated simulations. That is, we have 3 rates $\boldsymbol{\lambda}=[\lambda_1,\lambda_2,\lambda_3]$ each of duration 1 hour. To gather better statistics, I instead simulate a process with $\boldsymbol{\hat{\lambda}}=[\boldsymbol{\lambda},\boldsymbol{\lambda},\ldots,\boldsymbol{\lambda}]$, repeatedly cycling through the $\boldsymbol{\lambda}$ vector to allow larger sample sizes for each of the $\lambda$'s (while preserving their sequence and relative durations).
The first step produces a distribution of waiting times (in $\Lambda$ space) like this
which compares very well to the expected theoretical CDF (i.e. a unit-rate exponential).
The second step produces the final ($t$ space) arrival times. The empirical distribution of waiting times matches well with theory here as well:
This figure is a little bit more busy, but can be interpreted as follows. First, the upper panel shows the component CDFs associated with a homogeneous Poisson process at each of the three $\lambda$'s. We expect the aggregate CDF to be a mixture of these components
$$\Pr\big[\tau<\Delta{t}\big]=\sum_iw_i\Pr\big[\tau<\Delta{t}\mid\lambda=\lambda_i\big]$$
where $w_i=\Pr\big[\lambda=\lambda_i\big]$ are the mixing fractions.
For each component process, the expected number of samples is $\langle{n_i}\rangle=\lambda_iT_i$. Since the durations $T_i$ are equal, the expected mixing weights will then scale with the rates, i.e.
$$w_i\propto\lambda_i$$
The lower panel above shows the expected CDF obtained by mixing the components in proportion to their rates. As can be seen, the empirical CDF of the inhomogeneous simulation is consistent with this theoretically expected CDF.
Example Code
The following code can be run here. (Note: the algorithm is only 4 lines; the bulk of the code is comments and verification.)
%% SETUP
lam0=[1.5,2.1,3.4]; dt0=ones(size(lam0)); % rates and their durations
Nrep=1e3; lam=repmat(lam0,[1,Nrep]); dt=repmat(dt0,[1,Nrep]); % do replications (for stats check)
%% SIMULATION
L=cumsum([0,lam.*dt]); t=cumsum([0,dt]); % cumulative expected # events and time
Lmax=L(end); N=poissrnd(Lmax); % sample total # events
Lsmp=Lmax*sort(rand([1,N])); % sample event times (in "L space")
tsmp=interp1(L,t,Lsmp); % transform to "t space"
%% STATS CHECK
% "L space" waiting time CDF
dL=sort(diff(Lsmp)); p=(1:N-1)/N; % simulated
p0_L=1-exp(-dL); % exponential
h=plot(dL,p,'k',dL,p0_L,'r','LineWidth',1.5); set(h(1),'LineWidth',4);
axis tight; xlabel('dL'); legend('data','theory (exp)',4); title('L space CDF');
% "t space" waiting time CDF
dT=sort(diff(tsmp)); % simulated
wcmp=(lam0.*dt0)/(lam0*dt0'); pcmp=1-exp(-lam0'*dT); p0=wcmp*pcmp; % mixture
subplot(211); plot(dT,pcmp); ylabel('CDF'); title('Mixture Components'); axis tight;
legend(cellstr(num2str(lam0','lam = %0.1f ')),4);
subplot(212); h=plot(dT,p,'k',dT,p0,'r','LineWidth',1.5); set(h(1),'LineWidth',4);
axis tight; xlabel('dt'); ylabel('CDF'); title('Aggregate');
legend('data','theory (exp mixture)',4);
% mean arrival rate
tbin=0.5:2.5; Navg=hist(mod(tsmp,3),tbin)/Nrep; bar(tbin,Navg,1);
hold on; plot(tbin,lam0,'r.'); hold off;
legend('observed (hist.)','theory (lam.)',0);
xlabel('hour'); ylabel('arrivals/hour');
As requested in the comments, here is an example of the simulated mean arrival rates (see also updated code above):
Best Answer
Let $t=T_F$. Conditional on the number of occurences $N=n$, the arrival times $t_1,t_2,\dots,t_N$ are known to have the same distribution as the order statstics of $n$ iid unif$(0,t)$ random variables. Hence, the likelihood becomes \begin{align} L(\lambda,t) &= P(N=n) f(t_1,t_2,\dots,t_N|N=n) \\ &= \frac{e^{-\lambda t}(\lambda t)^n}{n!}\frac{n!}{t^n} \\ &= e^{-\lambda t}\lambda^n. \end{align} for $t\ge t_n$ and zero elsewhere. This is maximised for $\hat t=t_n$ and $\hat\lambda=n/t_n$. These MLEs don't exist if there are no occurrences $N=0$, however. Conditional on $N=n$, again using the fact that $t_n$ can be viewed as an order statistic (the maximum) of $n$ iid unif$(0,t)$ random variables, $E(t_N|N=n)=\frac n{n+1} t$. Hence, the estimator $t^*=\frac {n+1}n t_n$ is unbiased for $t$ conditional on $N=n$ and hence also conditional on $N\ge 1$. A reasonable frequentist estimator of $\lambda$ might be $\lambda^* = n/t^* = \frac{n^2}{(n+1)t_n}$ but this does not have finite expectation when $N=1$ so assessing its bias is even more troublesome.
Bayesian inference using independent, non-informative scale priors on $\lambda$ and $t$ on the other hand leads to a posterior $$ f(\lambda,t|t_1,\dots,t_N) \propto e^{-\lambda t}\lambda^{n-1}t^{-1}. $$ for $t>t_n,\lambda>0$. Integrating out $\lambda$, the marginal posterior of $t$ becomes $$ f(t|t_1,\dots,t_N) = \frac{n t_n^n}{t^{n+1}}, t>t_n, $$ and the posterior mean $E(t|t_1,\dots,t_N)=\frac n{n-1} t_n$. A $(1-\alpha)$-credible interval for $t$ is given by $\left(\frac{t_n}{(1-\alpha/2)^{1/n}}, \frac{t_n}{(\alpha/2)^{1/n}}\right)$.
The marginal posterior of $\lambda$, \begin{align} f(\lambda|t_1,\dots,t_N) &\propto \int_{t_\text{max}}^\infty e^{-\lambda t}\lambda^{n-1}t^{-1} dt \\ &= \lambda^{n-1}\Gamma(0,\lambda t_n) \end{align} where $\Gamma$ is the incomplete gamma function.