The Prothero-Robinson stiff ODE model

Bartel & Günther present a model for stiff systems as adapted from Kvaerno. The Prothero Robinson model that they use is defined by the following set of non-homogeneous equations:

 

The above equations represent sinusoidal solutions governed by the frequency w, with the stiffness factor governed by gamma. Bartel & Günther suggest -10E4 for the stiffness factor, with zero initial conditions at t=0, w=100 and eps =0.1. The solutions of the ODE’s are thus two variables with vastly differing settling rates, clearly justifying the classification as stiff. 

Plotted below is the solution obtained using the Generalized Rosenbrock Algorithm with the VELDS parameters.