How hot is the water in the pot? Question: How hot is the water in the pot? More precisely speaking, how can I get a temperature of the water as a function of time a priori? 
Background & My attempt: Recently I started spend some time on cooking. And I'm curious about it. I have learned mathematics as a undergraduate student for four years, but I know a little about thermodynamics. (I listened to such a lecture once. So I've heard of $dU = TdS - pdV$, Entropy and Gibbs energy for example though I forgot almost everything; anyway I think I've never seen a formula depending on time.)   So I conduct a small experiment first: I heat 100ml of water by IH correspond approximately to 700W and measure its temperature every 30 seconds. Here is the results.

It looks almost linear, but I think linear approximation is inappropriate; because if so, the water gets higher than $100^\circ\mathrm{C}$. So I guess it's some convex increasing function like $T(t) = 100 - \alpha e^{-t/\beta}$ for some positive constant $\alpha$ and $\beta$. But it doesn't fit the data.  (It does fit the data. I just made a mistake in simple calculation. See my answer.)
I think I ignored too many factors. So feel free to assume anything reasonable. I would greatly appreciate if you help me. Thank you.
Additional question: I do a experiment and some calculation to deal with a problem pointed out in the comments of my answer: bad fitting at lower temperature. However, I cannot get a better solution. Fitting seems worse than before... Here is the results what I got. I heated 100ml water in pot with 9cm radius by IH correspond to 700W. (For calculation, I added linear interpolation values in the graph.) How can I get a better solution?

(Light blue curve is a logistic approximation defined by $T = \dfrac{100}{1 + 1.62 e^{-0.0168 t}}$ as mentioned here.)
 A: It seems we have reached the point where simple models are no longer satisfying. Rather than posing ad hoc DEs maybe it's time to try an actual physical model. Short of doing a full hydrodynamic simulation (definitely overkill here) we can try what is called a lumped capacitance model where we divide the system up into a number of "lumps" and energy flows between the lumps:

The fundamental law governing this system is the conservation of energy. Every lump has an equation of the form
$$ \frac{\mathrm{d}}{\mathrm{d}t}(\text{energy in lump}) = \text{rate of energy entering} - \text{rate of energy leaving}. $$
We treat the heat input as a fixed flow and the air (environment) as a heat bath at a fixed temperature. If we let the heat capacity of the $i$-th lump be $C_i(T)$, which can be a function of temperature then
$$ \begin{array}{rcl}
\frac{\mathrm{d}}{\mathrm{d}t} \left( C_p(T_p) T_p \right) &=& P - r_1 (T_p - T_w) - r_2 (T_p - T_a)\\
\frac{\mathrm{d}}{\mathrm{d}t} \left( C_w(T_w) T_w \right) &=& r_1 (T_p - T_w) - r_3 (T_w - T_a)
\end{array} $$
You can look up how the heat capacity $C_w(T)$ varies with temperature for water (though probably not for the pot material?), but we will simplify dramatically and assume, incorrectly, that the heat capacities are constant.
$$ \begin{array}{rccl}
\frac{\mathrm{d}}{\mathrm{d}t} T_p &=& - \frac{r_1 + r_2}{C_p} T_p + \frac{r_1}{C_p} T_w + \frac{P + r_2 T_a}{C_p} &\equiv&  a T_p + b T_w + s_p\\
\frac{\mathrm{d}}{\mathrm{d}t} T_w &=& \frac{r_1}{C_w} T_p - \frac{r_1 + r_3}{C_w} T_w + \frac{r_3 T_a}{C_w} &\equiv& c T_p + d T_w + s_w
\end{array} $$
where the $a,b,c,d,s_p,s_w$ are shorthands. Note that only six combinations of the seven parameters ($r_1,r_2,r_3,P,T_a,C_p,C_w$) actually enter the problem, so there is some degeneracy of the parameters. You can see, Taro, that this is almost the model you came up with in your answer. The difference is that I'm including the heat input explicitly, so that conservation of energy is guaranteed.
With the obvious matrix shorthand these equations can be written
$$ \dot{T}-MT = s, $$
which, for a constant source, has the solution
$$ T\left(t\right) = \mathrm{e}^{Mt}T_{0}+\mathrm{e}^{Mt}\left(\int_{0}^{t}\mathrm{e}^{-M\tau}\mathrm{d}\tau\right) s. $$
When $M$ is invertible (which it definitely should be for this problem - if it's not then there is a mistake somewhere) the integral can be simplified:
$$ T\left(t\right) = \mathrm{e}^{Mt}\left(T_{0}+M^{-1}s\right)-M^{-1}s. $$
You can check this satisfies the original equation with the proper boundary conditions. There are eight parameters to fit: the four matrix elements, two sources, and two initial temperatures. It is a non-linear regression problem since the matrix exponential depends non-linearly on the fit parameters. So I'm afraid I don't know of a robust and efficient way to fit this model to your data, unless you use some assumptions to simplify the parameter dependence, but this is the physically motivated model for your situation.
A: To close this post, I write an answer by myself though it turned out that I made just a simple calculation mistake.
Since the temperature increase should be monotonic and approach to zero at boiling point, it's reasonable to assume that the temperature increase $dT/dt$ is proportional to the difference $T - 100 \mathrm{^\circ C}$, that is, 
$$ \frac{dT}{dt} = -k(T - 100 \mathrm{^\circ C}) $$
holds for some positive constant $k$.
Solving this equation gives 
$$ T = 100 \mathrm{^\circ C} + (T(t_0) - 100 \mathrm{^\circ C})e^{-k(t - t_0)}. $$
Let's determine coefficient $k$ from $N$ measurements by linear regression. Let $c$ be a time interval of experiments and $x_n$ the temperature of the water on $t_n = cn$. Then estimate the slope of the tangent line by
$$ y_n = \mathrm{mean}\big(\frac{t_{n + 1} - t_{n}}{c}, \frac{t_{n} - t_{n - 1}}{c}\big) = \frac{t_{n + 1} - t_{n - 1}}{2c} $$
for $0 < n < N$. 
From above equation, there should be a relation of the form
$$ y_n = -k(x_n - 100 \mathrm{^\circ C}) + \varepsilon_n $$
where $\varepsilon_n$ stands for experimental errors. I denote this equation by
$$ y = -kx +\varepsilon $$
as shorthand. The best estimator $\hat{k}$ is given if $x$ and $\varepsilon$ are orthogonal to each other. Therefore 
$$ \hat{k} = -\frac{(x, y)}{(x, x)} \approx 0.00667 \mathrm{s^{-1}} $$
from calculations. And this value fits the experimental data well.

Note:
I completely rewrite this answer. Here, I would like to review where my answer was inappropriate. It seems not a problem of physic but a problem of statistics. Last time, I solved DE first and took logarithm to make it linear. However the experimental errors also transformed. Especially, $\ln(100 - x_n) \to -\infty$ as $x_n \to 100$. So this seems to cause overfitting at higher temperature and bad fitting at lower temperature. (Considering a effect of the pot looks a good idea but everything I tried fails. It won't fit data and still open though I already got a reasonable approximation.) 
Thank you so much for helping me, Chris, Stefan Bischof, Michael Brown and Christoph. 
A: Data fitting usually is done by the least-squares fit method. It's already implemented in a couple of computer algebra systems. 
I recommand the open source gnuplot for fitting the data. You will find the most appropriate model results in the smallest number of the least squares fit. My link contains an gnuplot example script for fitting.
Answer
$$T_1(t) = A\cdot t - B\cdot e^{-k\cdot t}+C$$
Start fitting with the suggestion of saturated increase of Chris in the comments. Try to identify $T_{\text{boil}}(p)\approx 100\,$°C and $T_{\text{room}}$ from $T(0s)$. 
Edit about extended time after your measurement:
Keeping your hot plate still on lets me think about the fact that the water is still cooking. Conversion from water to gas phase is still in process. A measurement diagram after the heating phase would look like a constant/noisy $\approx 100\,$°C saturation. This extended measurement better would be covered by some hyperbolic tangent 
$$T_2(t)=A\cdot t+ B\cdot \tanh(k\cdot t)+C$$
The logistic growth represents slower temperature increase and asymptotical increase to $\approx 100\,$°C.
