% An Example of stiff ODE system:
%     y1' = -49.5*y1 + 50.5*y2
%     y2' =  50.5*y1 - 49.5*y2
% with initial conditions y1(0)=1, y2(0)=-1.
% We need to obtain an approximate solution for 0<=t<=T.
% Try to enter different values of N to see what happens.
% First try relatively large N, say N>=10. You will see that both the
% explicit and implicit methods generate similar approximations. When
% you reduce N to 4, the difference is shocking. The explicit method
% got it totally wrong, while the implicit method maintains a decent
% approximation.
T = 0.1;
N = input('enter number of subintervals N: ');
tgrid = linspace(0,T,N+1);
h = T/N;
w = zeros(2,N+1);
w(1,1) = 1; w(2,1)=-1;
A = [-49.5 50.5; 50.5 -49.5];
lambda = max(abs(eig(A)));
fprintf('explicit factor is %6.4f\n',1-lambda*h)
B = eye(2) - h*A;
% First we use the explicit Euler:
for j=1:N
    w(:,j+1) = w(:,j) + h*A*w(:,j);
end
figure
plot(tgrid,w(1,:),'-',tgrid,w(2,:),'--');
title('Solution from explicit Euler');
legend('y1','y2');

% Then we use the implicit Euler:
for j=1:N
    w(:,j+1) = B\w(:,j);
end
figure
plot(tgrid,w(1,:),'-',tgrid,w(2,:),'--');
title('Solution from implicit Euler');
legend('y1','y2');
