% Problem specification.
M = [-10 -10; 40 -10];
g = @(t) [sin(t); cos(t)];
f = @(t, y) M * y + g(t);
tint = [0, 10];
y0 = [5210; 20259] / 249401;

% Simple implementation of the backward Euler method.
N = 100; % Vary this and see that there is no instability.
h = diff(tint) / N;

t = linspace(tint(1), tint(2), N+1);
y = zeros(2, N+1);
y(:, 1) = y0;

invmat = inv(eye(2) - h*M);
for n = 1:N
   y(:, n+1) = invmat * (y(:, n) + h * g(t(n+1)));
end

plot(t, y)
