% Explicit RK4 Stiff
%
% u'' + w*u = 0
%
w=1; mu=-100
dt=0.02; tmax = 30;
n=tmax/dt;
A = [0 w 0; -w 0 0; 0 1 mu];
Id = eye(3);
u = zeros(3,n+1); u(:,1)= [1; 0; 0];
B=(A*dt+Id)
for j=1:n
    u(:,j+1) = B*u(:,j);
end
t=[0:dt:tmax];
plot(t,u(1,:),t,u(3,:),t,cos(t)); xlabel('t'); ylabel('u');