% Implicit Euler Linear Pendulum
%
% u'' + w*u = 0
%
w=1; mu=-100
dt=0.03; 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];
B1=0.5*A*dt+Id;
B2=A*dt+Id;
for j=1:n
    s1 = A*u(:,j)*dt;
    s2 = A*(u(:,j)+0.5*s1)*dt;
    s3 = A*(u(:,j)+0.5*s2)*dt;
    s4 = A*(u(:,j)+s3)*dt;
    u(:,j+1) = u(:,j)+(s1+2*s2+2*s3+s4)/6;
end
t=[0:dt:tmax];
plot(t,u(1,:),t,u(3,:),t,cos(t)); xlabel('t'); ylabel('u');