Settling time of an op-amp integrator
First, I will present a method that uses Mathematica to solve this problem. When I was studying this stuff I used the method all the time (without using Mathematica of course). Besides that, this answer will be using 'heavy' mathematics ;).

Well, we are trying to analyze the following opamp-circuit:



simulate this circuit &ndash; Schematic created using CircuitLab
When we use and apply KCL, we can write the following set of equations:
$$
\begin{cases}
\text{I}_1=\text{I}_\text{k}+\text{I}_5\\
\\
\text{I}_5=\text{I}_2+\text{I}_3\\
\\
0=\text{I}_2+\text{I}_3+\text{I}_6\\
\\
\text{I}_\text{o}+\text{I}_\text{p}+\text{I}_1=\text{I}_4
\end{cases}\tag1
$$
When we use and apply Ohm's law, we can write the following set of equations:
$$
\begin{cases}
\text{I}_1=\frac{\text{V}_\text{i}-\text{V}_1}{\text{R}_1}\\
\\
\text{I}_2=\frac{\text{V}_1-\text{V}_3}{\text{R}_2}\\
\\
\text{I}_3=\frac{\text{V}_1-\text{V}_3}{\text{R}_3}\\
\\
\text{I}_4=\frac{\text{V}_3}{\text{R}_4}\\
\\
\text{I}_\text{p}=\frac{\text{V}_2}{\text{R}_5}
\end{cases}\tag2
$$
Substitute \$(2)\$ into \$(1)\$, in order to get:
$$
\begin{cases}
\frac{\text{V}_\text{i}-\text{V}_1}{\text{R}_1}=\text{I}_\text{k}+\frac{\text{V}_2}{\text{R}_5}\\
\\
\frac{\text{V}_2}{\text{R}_5}=\frac{\text{V}_1-\text{V}_3}{\text{R}_2}+\frac{\text{V}_1-\text{V}_3}{\text{R}_3}\\
\\
0=\frac{\text{V}_1-\text{V}_3}{\text{R}_2}+\frac{\text{V}_1-\text{V}_3}{\text{R}_3}+\text{I}_6\\
\\
\text{I}_\text{o}+\text{I}_\text{p}+\frac{\text{V}_\text{i}-\text{V}_1}{\text{R}_1}=\frac{\text{V}_3}{\text{R}_4}
\end{cases}\tag3
$$
Now, using an ideal opamp, we know that:

$$\text{V}_x:=\text{V}_+=\text{V}_-=\text{V}_1=\text{V}_2\tag4$$
$$\text{I}_\text{p}=\text{I}_\text{k}=0\space\text{A}\tag5$$

So we can rewrite equation \$(3)\$ as follows:
$$
\begin{cases}
\frac{\text{V}_\text{i}-\text{V}_x}{\text{R}_1}=0+\frac{\text{V}_x}{\text{R}_5}\\
\\
\frac{\text{V}_x}{\text{R}_5}=\frac{\text{V}_x-\text{V}_3}{\text{R}_2}+\frac{\text{V}_x-\text{V}_3}{\text{R}_3}\\
\\
0=\frac{\text{V}_x-\text{V}_3}{\text{R}_2}+\frac{\text{V}_x-\text{V}_3}{\text{R}_3}+\text{I}_6\\
\\
\text{I}_\text{o}+0+\frac{\text{V}_\text{i}-\text{V}_x}{\text{R}_1}=\frac{\text{V}_3}{\text{R}_4}
\end{cases}\tag6
$$
Now, we can solve for the transfer function:
$$\mathscr{H}:=\frac{\text{V}_3}{\text{V}_\text{i}}=-\frac{\text{R}_2\text{R}_3}{\text{R}_1\left(\text{R}_2+\text{R}_3\right)}\tag7$$
Where I used the following Mathematica-code:
In[1]:=V1 = Vx;
V2 = Vx;
Ik = 0;
Ip = 0;
FullSimplify[
 Solve[{I1 == Ik + I5, I5 == I2 + I3, 0 == I2 + I3 + I6, 
   Io + Ip + I1 == I4, I1 == (Vi - V1)/R1, I2 == (V1 - V3)/R2, 
   I3 == (V1 - V3)/R3, I4 == V3/R4, Ip == V2/R5}, {Io, I1, I2, I3, I4,
    I5, I6, Vx, V3}]]

Out[1]={{Io -&gt; -(((R3 R4 + R2 (R3 + R4)) Vi)/(R1 (R2 + R3) R4)), I1 -&gt; Vi/R1,
   I2 -&gt; (R3 Vi)/(R1 (R2 + R3)), I3 -&gt; (R2 Vi)/(R1 (R2 + R3)), 
  I4 -&gt; -((R2 R3 Vi)/(R1 R2 R4 + R1 R3 R4)), I5 -&gt; Vi/R1, 
  I6 -&gt; -(Vi/R1), Vx -&gt; 0, V3 -&gt; -((R2 R3 Vi)/(R1 (R2 + R3)))}}

My equation was also confirmed using LTspice.

When we want to apply the derivation from above to your circuit we need to use Laplace transform (I will use lower case function names for the functions that are in the (complex) s-domain, so \$\text{y}\left(\text{s}\right)\$ is the Laplace transform of the function \$\text{Y}\left(t\right)\$):
$$\text{R}_2=\frac{1}{\text{sC}}\tag8$$
So, we can rewrite the transfer function as:
$$\mathcal{h}\left(\text{s}\right)=-\frac{\frac{1}{\text{sC}}\cdot\text{R}_3}{\text{R}_1\left(\frac{1}{\text{sC}}+\text{R}_3\right)}=-\frac{\text{R}_3}{\text{R}_1\left(1+\text{sCR}_3\right)}=-\frac{\text{R}_3}{\text{R}_1}\cdot\frac{1}{1+\text{sCR}_3}\tag9$$
Now, the function of a square wave is given by:
$$\text{V}_\text{i}\left(t\right)=\frac{4\hat{\text{u}}}{\pi}\sum_{\text{n}\ge1}\frac{\sin\left(2\pi\left(2\text{n}-1\right)\text{f}t\right)}{2\text{n}-1}\tag{10}$$
Where \$\hat{\text{u}}\$ is the amplitude of the square wave and \$\text{f}\$ is the frequency of the square wave.
Using the convolution property of the Laplace transform, we can write:
$$\text{V}_3\left(t\right)=\int_0^t\mathscr{L}_\text{s}^{-1}\left[\text{v}_\text{i}\left(\text{s}\right)\right]_{\left(t-\tau\right)}\cdot\mathscr{L}_\text{s}^{-1}\left[-\frac{\text{R}_3}{\text{R}_1}\cdot\frac{1}{1+\text{sCR}_3}\right]_{\left(\tau\right)}\space\text{d}\tau=$$
$$-\frac{\text{R}_3}{\text{R}_1}\int_0^t\text{V}_\text{i}\left(t-\tau\right)\cdot\mathscr{L}_\text{s}^{-1}\left[\frac{1}{1+\text{sCR}_3}\right]_{\left(\tau\right)}\space\text{d}\tau\tag{11}$$
Using the table of selected Laplace transforms, we can write:
$$\text{V}_3\left(t\right)=-\frac{\text{R}_3}{\text{R}_1}\int_0^t\text{V}_\text{i}\left(t-\tau\right)\cdot\frac{\exp\left(-\frac{\tau}{\text{CR}_3}\right)}{\text{CR}_3}\space\text{d}\tau=$$
$$-\frac{\text{R}_3}{\text{R}_1}\cdot\frac{1}{\text{CR}_3}\int_0^t\text{V}_\text{i}\left(t-\tau\right)\cdot\exp\left(-\frac{\tau}{\text{CR}_3}\right)\space\text{d}\tau=$$
$$-\frac{1}{\text{CR}_1}\int_0^t\text{V}_\text{i}\left(t-\tau\right)\cdot\exp\left(-\frac{\tau}{\text{CR}_3}\right)\space\text{d}\tau\tag{12}$$
Using \$(10)\$ we can write:
$$\text{V}_3\left(t\right)=-\frac{1}{\text{CR}_1}\int_0^t\left\{\frac{4\hat{\text{u}}}{\pi}\sum_{\text{n}\ge1}\frac{\sin\left(2\pi\left(2\text{n}-1\right)\left(t-\tau\right)\text{f}\right)}{2\text{n}-1}\right\}\cdot\exp\left(-\frac{\tau}{\text{CR}_3}\right)\space\text{d}\tau=$$
$$-\frac{4\hat{\text{u}}}{\text{CR}_1\pi}\sum_{\text{n}\ge1}\frac{1}{2\text{n}-1}\int_0^t\sin\left(2\pi\left(2\text{n}-1\right)\left(t-\tau\right)\text{f}\right)\cdot\exp\left(-\frac{\tau}{\text{CR}_3}\right)\space\text{d}\tau=$$
$$\frac{4\hat{\text{u}}}{\text{CR}_1\pi}\sum_{\text{n}\ge1}\frac{1}{1-2\text{n}}\int_0^t\sin\left(2\pi\left(2\text{n}-1\right)\left(t-\tau\right)\text{f}\right)\cdot\exp\left(-\frac{\tau}{\text{CR}_3}\right)\space\text{d}\tau\tag{13}$$

For now I am writing a code that gives your circuit, so that we can do some extra analysis on your circuit specific. That code is running, if I finish it I will publish it here in my answer.