%% Question 3 - Comparison of Implict & Explicit Euler

%% Initial setup
opt = odeset('RelTol',1e-5,'AbsTol',1e-6); % increase precision
t0 = 0;
T = 0.1;
x0 = [2; 1];
[t,x] = ode23(@linsystem, [t0,T], x0, opt);

%% Timestep $\tau = 0.002$
figure;
plot3(t, x(:,1), x(:,2), 'k', 'DisplayName', 'ode23');
hold on;
xlabel('t');
ylabel('x_1');
zlabel('x_2');
view(-127.5, 30);
xlim([0, 0.1]);
ylim([0 10]);
zlim([-8 2]);
h=0.002;
[te,xe] = eul(@linsystem, t0, T, x0, h);
plot3(te, xe(:,1), xe(:,2), '-r', 'DisplayName', 'Euler');
[ti,xi] = ieuler(@linsystem, t0, T, x0, h);
plot3(ti, xi(:,1), xi(:,2), '-b', 'DisplayName', 'Implicit Euler (Fixed Point)');
[tin,xin] = ieuler_newton(@linsystem_newton, t0, T, x0, h);
plot3(tin, xin(:,1), xin(:,2), '-m', 'DisplayName', 'Implicit Euler (Newton)');
legend('Location', 'NorthWest');

%% Timestep $\tau = 0.0021$
figure;
plot3(t, x(:,1), x(:,2), 'k', 'DisplayName', 'ode23');
hold on;
xlabel('t');
ylabel('x_1');
zlabel('x_2');
view(-127.5, 30);
xlim([0, 0.1]);
ylim([0 10]);
zlim([-8 2]);
h=0.0021;
[te,xe] = eul(@linsystem, t0, T, x0, h);
plot3(te, xe(:,1), xe(:,2), '-r', 'DisplayName', 'Euler');
[ti,xi] = ieuler(@linsystem, t0, T, x0, h);
plot3(ti, xi(:,1), xi(:,2), '-b', 'DisplayName', 'Implicit Euler (Fixed Point)');
[tin,xin] = ieuler_newton(@linsystem_newton, t0, T, x0, h);
plot3(tin, xin(:,1), xin(:,2), '-m', 'DisplayName', 'Implicit Euler (Newton)');
legend('Location', 'NorthWest');

%% Timestep $\tau = 0.0019$
figure;
plot3(t, x(:,1), x(:,2), 'k', 'DisplayName', 'ode23');
hold on;
xlabel('t');
ylabel('x_1');
zlabel('x_2');
view(-127.5, 30);
xlim([0, 0.1]);
ylim([0 10]);
zlim([-8 2]);
h=0.0019;
[te,xe] = eul(@linsystem, t0, T, x0, h);
plot3(te, xe(:,1), xe(:,2), '-r', 'DisplayName', 'Euler');
[ti,xi] = ieuler(@linsystem, t0, T, x0, h);
plot3(ti, xi(:,1), xi(:,2), '-b', 'DisplayName', 'Implicit Euler (Fixed Point)');
[tin,xin] = ieuler_newton(@linsystem_newton, t0, T, x0, h);
plot3(tin, xin(:,1), xin(:,2), '-m', 'DisplayName', 'Implicit Euler (Newton)');
legend('Location', 'NorthWest');

%% Timestep $\tau = 0.0002$
figure;
plot3(t, x(:,1), x(:,2), 'k', 'DisplayName', 'ode23');
hold on;
xlabel('t');
ylabel('x_1');
zlabel('x_2');
view(-127.5, 30);
xlim([0, 0.1]);
ylim([0 10]);
zlim([-8 2]);
h=0.0002;
[te,xe] = eul(@linsystem, t0, T, x0, h);
plot3(te, xe(:,1), xe(:,2), '-r', 'DisplayName', 'Euler');
[ti,xi] = ieuler(@linsystem, t0, T, x0, h);
plot3(ti, xi(:,1), xi(:,2), '-b', 'DisplayName', 'Implicit Euler (Fixed Point)');
[tin,xin] = ieuler_newton(@linsystem_newton, t0, T, x0, h);
plot3(tin, xin(:,1), xin(:,2), '-m', 'DisplayName', 'Implicit Euler (Newton)');
legend('Location', 'NorthWest');
