function fn = oscil(t, y)
    m = 1; c = 2; k = 25;
    omega = 1;
%    omega = sqrt(k/m - (c/m)^2/2);
    [n,j] = size(y); fn = zeros(n,j);
    fn(1) = y(2);
    fn(2) = -c*y(2)/m - k*y(1)/m + (20*cos(omega*t) + 52*sin(omega*t))/m;
return