Band pass filter unstable and lower gain than simulation

In my first answer I have described how you can find the exact solution for the 2 zeros and the 2 pole frequencies (which are identical to the wanted break frequencies). However, here is a good approach which should be sufficient for the shown circuit. In principle, I follow the way as outlined already by Dave Tweed´s answer: Simplification of the circuit. In the present case, you can create three different (simplified) circuits of first order only which easily can be analysed.

1.) For the first rising region of the transfer function the high pass part with C1 is responsible (C2 causes the falling part and can be neglected). Furthermore, for very low frequencies (including DC) the gain Ao=1+R3/R2 is assumed to be not much larger than unity which is the possible minimum.

Hence, for acceptable filtering it is assumed that R2>>R3. As an equivalent diagram for the lower frequency range (without C2 and R2) we arrive at a circuit with only the three components R1, R3 and C1. It is a simple task to find the relevant time constants (invers to the corresponding break frequencies):
Using your indices, we thus find T2=(R1+R3)C1 and T1=R1C1.

2.) Above the frequency f1 the capacitor C1 is not effective any more (and the capacitor C2 is assumed to be not yet effective). Hence, we have a simple non-inverting amplifier with the gain (maximum of the transfer function) Amax= 1+R3/Rp with Rp=R1||R2.

3.) For rising frequencies, the low pass part with capacitor C2 becomes effective (C1 is considered as a short). Hence, the feedback path consists of R3||C2 and Rp only.

The time constant T3 (pole frequency) can be derived as T3=R3C2 and the last break frequency (zero) is determined by T4=R3C2/(1+R3/Rp).

Finally, it is to be noted, that all results are in agreement with the values given in the scetched BODE diagram. This can be verified using the well-known relations for a 20dB gain slope (as used in the graph with G1/G2=f1/f2).

Final remark: Thus, it can be concluded that the information contained in the scetched BODE diagram (break frequencies) also are only approximations.

That's the problem with simulating with "ideal" components — you see behaviors that you would never see in the real world.

Your circuit has no resistance in it anywhere. The function generator is an ideal voltage source with zero output resistance. The oscilloscope has infinite input resistance (open circuit). And the components have no parasitic series or parallel resistance, either.

Therefore, the behavior that the simulator is showing you is correct. The 503 Hz sine wave is the L-C circuit continuing to "ring" from the start-up transient. This ringing will never die out. And you see none of the 5 kHz square wave at the output because your filter has infinite "Q" (quality factor), which means it blocks other frequencies perfectly.

When simulating a circuit that has only ideal components, you need to remember to model the parasitic effects of real components. Depending on the accuracy that you need, you might include the series resistance of both inductors and capacitors, and maybe some parallel capacitance on inductors and parallel resistance on capacitors as well. Usually, when simulating more complex circuits — which almost always contain resistors anyway — the effects of these parasitic components are insignificant.