The first thing you should have done is to put your schematic into the schematic editor here. It's a lot easier on the eye:
simulate this circuit – Schematic created using CircuitLab
(I held to your use of \$U\$ instead of \$V\$.)
I'd use brute-force KCL in the following way (using Python in SageMath and SymPy):
eq1 = Eq( u2, u3 ) # Ideal opamp
eq2 = Eq( u2*s*c1 + u2/r1, u1*s*c1 + u4/r1 ) # KCL U2
eq3 = Eq( u3*s*c2, io + u4*s*c2 ) # KCL U3
eq4 = Eq( u4/r1 + u4/r2 + u4*s*c2, u2/r1 + u3*s*c2 ) # KCL U4
ans = solve( [eq1, eq2, eq3, eq4], [u2, u3, u4, io] ) # solve system
i1 = simplify( (ans[u2]-u1)*s*c1 )
zin = simplify( u1/i1 )
At this point, I do have the result. You can substitute \$s=j\omega\$ into the result, if needed. And/or try out some real part values, as well, and run a simulation to validate the approach. I'd recommend doing that. For example:
simulate this circuit
Theory says:
abs( zin.subs( {s:I*2*pi*100, r1:5e3, r2:2e3, c1:1e-6, c2:3e-6} ) ).n()
18623.6083515058
Placing this circuit into LTspice, using an LT1800 opamp model, says:
So it shows only a tiny difference out at the sixth decimal place. I'd call that "close enough" for validation of the approach.
Added:
zin
-c2*r1*r2*s - r1 - r2 - 1/(c1*s)
abs(zin.subs(s,I*omega))
sqrt(c2**2*omega**2*r1**2*r2**2 + r1**2 + 2*r1*r2 + r2**2 - 2*c2*r1*r2/c1 + 1/(c1**2*omega**2))
