A stronger version of Cauchy-Goursat Theorem Claim: $C$ is a simple closed contour, $f$ is continuous at all points interior and on $C$, and $f$ is analytic at all points interior to $C$, then $\int_C f(z)dz = 0$.
To prove this, I suppose we may somehow approach $\int _{C}f(z) dz$ by a sequence contour integrals on closed curve inside $C$, which has value zero by the usual Cauchy Goursat theorem, but can't quite finish the proof.
I saw a similar quesiton Here, but not quite satisfied with the answer, in particular, the step of bounding the difference of integrals around two different contours. I think the post refers to the ML theorem to bound each contour integrals, but I am not sure how the separate bound for each can be put together to give a bound for the difference.
And ideas or comments are much appreciated.
 A: *

*If $C$ is supposed rectifiable, your claim is correct, but the proof is non-trivial. See for example "The Cauchy-Goursat theorem for rectifiable jordan curves" by Walsh. However, if you know Mergeylan's  theorem (for example, Rudin, in "Real&Complex analysis" provides a proof of it in the last chapter), the result is quite easy to obtain: the claim is trivially true for any polynomial function, and thanks to Mergeylan's theorem, every holomorphic function on $\text{Int}(C)$ continuous on $C$ can be uniformly approximated by a sequence of polynomial on $\text{Int}(C)\cup C$, and the result follows.

*If $C$ is not supposed to be rectifiable, your claim is wrong, as the integral $\int_C f(z)$ could even be undefined (an example is presented here)

*If $C$ is supposed to be rectifiable convex curve, there is an elementary proof sketched here, althought it is quite unelegant.


Counterexample with a non-rectifiable curve:
$$\gamma:[0,1]\to \mathbb{R}^2\\
\gamma(t)=\left(t,t\sin\left(\frac{1}{t^2}\right)\right)$$
This contour is simple and continuous (adding a path connecting the ends we can
construct also a closed contour). However, $$\int_\gamma zdz$$ does not exist as a limit of riemann sums, as one can show without much effort (hint: consider two suitable Riemann sums and prove that their difference does not go to $0$).
A more pathological counterexample is obtained by choosing as contour a closed osgood curve, which has positive area.
Proof if $C$ is rectifiable and convex
If $C$ is rectifiabe, let us call $C_n$ a polygonal path that approximates the curve (and has length $l\le \text{length}(C)$), defined as follows: let $\gamma:[0,1]\to C$ be a parametrization of the contour. Given a partition $\mathcal{P}$of $[0,1]$, the associated polygonal is defined as the polygonal line connecting $\gamma(0),\gamma(t_1),\dots,\gamma(t_{n-1}),\gamma(1)$.
In the following, we will write $\gamma_{[t_i,t_{i+1}]}$ for the curve with domain restricted to $[t_i,t_{i+1}]$, $\gamma(t_i)\to \gamma(t_{i+1})$ to represent the straight segment connecting the two points, and $\Delta_i$ for the length of such segment. We will also associate to $\mathcal{P}$ a Riemann sum $\mathcal{S}=\sum f(\gamma(t_{*}))\cdot (\gamma(t_{k+1})-\gamma(t_k))$ (where $t_*\in [t_i,t_{i+1}]$) Now let us analyze
$$\left|\int_C f(z)dz-\int_{C_n}f(z)dz\right|$$
By triangle inequality, we have
\begin{align*}\left|\int_C f(z)dz-\int_{C_n}f(z)dz\right|&\le\left|\int_C f(z)dz-\mathcal{S}\right|+
\left|\mathcal{S}-\int_{C_n}f(z)dz\right|\\
\left|\mathcal{S}-\int_{C_n}f(z)dz\right|&\le\sum\int_{\gamma(t_i)\to \gamma(t_{i+1})}|f(z)-f(t_*)|dz\\&\le \sum \left(\max_{\gamma(t_i)\to \gamma(t_{i+1})}(f)-
\min_{\gamma(t_i)\to \gamma(t_{i+1})}(f)\right)\Delta_i
\end{align*}
Since $f$ is continuous on $C$, which is bounded and closed, $f$ is uniformly continuous on $C$. Since it is holomorphic in the interior, it is uniformly continuous on $C_n$. Also, since $C$ is rectifiable, we know that $\int_\gamma f-S\to 0$ (since a continuous function on a rectifiable curve is integrable). Thus we can choose a partition with mesh small enough such that:
$$
\left|\int_C f(z)dz-\mathcal{S}\right|\le \frac{\varepsilon}{2}\\
\max_{\gamma(t_i)\to \gamma(t_{i+1})}(f)-
\min_{\gamma(t_i)\to \gamma(t_{i+1})}(f)\le \frac{\varepsilon}{2\text{length}(\gamma)}
$$
Thus
$$\left|\int_C f(z)dz-\int_{C_n}f(z)dz\right|\le \varepsilon$$
Since $\varepsilon$ was an arbitrary positive value, their difference must be $0$, i.e.
$$\int_C f(z)dz=\int_{C_n}f(z)dz$$
It remains only to prove that $\int_{C_n}f(z)dz=0$, which is easy: we will prove it by induction on $n$.
$n=3$: left as exercise
$n+1$: we can decompose $\oint C_n$ as
$$\int_{C_n}f(z)dz=\int_{\gamma(0)\to \gamma(t_1)\to \gamma(t_2)\to \gamma(0)}f(z)dz+
\int_{\gamma(t_2)\to \dots\to \gamma(1)\to \gamma(t_2)}f(z)dz$$
The first term on the RHS is zero by the base step, the second by induction hypotesis (since it is a polygonal with $n-1$ sides.
