Looks pretty good. I suspect you're getting the diverging solution because the Q term in the denominator of your resonator TF is not imaginary. This means that the resonator has no damping, and any drive would indeed lead to diverging amplitude.

You might need to define the TF from the command line function, rather than in the simulink block; I think the block throws an error when you give it complex coefficients. This is surprising to me, so maybe I'm missing something.

Edit: That's not it, I just wasn't thinking in s-domain.

Here are a few other things I notice about this model that might be important:

• The first BP filter is specified in Hz, rather than rad/s. This means you're inadvertently BPing well away from the plant's resonance.
• You're computing the 'amplitude' by low passing the bandpassed signal. But what's done in the model is an RMS followed by lowpass. That is, the squared signal is scaled and lowpassed. Otherwise, you should get zero (BP near 1kHz followed by LP at 10s of Hz, with no nonlinear steps => 0)
• You don't have the loop shaping filter, typically included in the PWRCTRL filter module. This contains the integrator needed to reach the correct setpoint (on long time scales, else you'll reach some offset from the setpoint) and an extra zero to increase the UGF.
• The butterworth LP on the lower path has a pretty high cutoff, usually I set this at a few to 10 Hz.
• The saturation gives a clean square wave iff there is a large gain just before, otherwise it gives a square wave with some rise/fall time where the original sine was crossing zero.
  • This is the source of at least some problems. Whenever I was increasing this gain, I found the output of differentiator decreased proportionally. I was also eventually getting some errors from the saturation block about detecting too many zero crossings. Turning off zero crossing detection at the saturation block solved both of these issues. In retrospect, this makes sense -- zero crossing detection is meant to smooth out discontinuities by decreasing the step size until they go away; this block is intentionally introducing such a discontinuity.
• The laplace transform of an exponentially decaying sine has all positive real coefficients in the denominator. While the frequency-domain transfer function is what Rana posted above and has a (omega_0^2-omega^2) term, in the laplace domain it's different. I'd guess this is why you're geting a divergence.