%% Chapter 4, Examples 1 & 2
% Explicit and Implicit Euler

%% Cleanup
close all;
clear f t_true y_true y1 h t1 i y2 t2 y3 t3 y4 t4;

%% Analytic solution
f = @(x) exp(-x/2);
t_true = linspace(0,20,50);
y_true = f(t_true);

%% Numerical solution for h=1
y1(1) = 1;
h = 1;
t1(1) = 0;
i = 1;
while (t1(i) <= 30)
    y1(i+1) = y1(i) - .5*h*y1(i);
    t1(i+1) = t1(i) + h;
    i = i + 1;
end

%% Numerical solution for h=4.2
y2(1) = 1;
h = 4.2;
t2(1) = 0;
i = 1;

while (t2(i) <= 30)
    y2(i+1) = y2(i) - .5*h*y2(i);
    t2(i+1) = t2(i) + h;
    i = i + 1;
end

%% Plot explicit solutions
figure(1); clf;
plot(t_true,y_true)
hold on
plot(t1,y1,'k.--','MarkerSize',20)
plot(t2,y2,'ro--','MarkerSize',10)
axis([0 20 -1.5 2])
xlabel('t')
ylabel('y(t)')
legend('Exact','Explicit Euler, h=1','Explicit Euler, h=4.2','Location','NorthWest')

%% Implicit solution for h=1
y3(1) = 1;
h = 1;
t3(1) = 0;
i = 1;
while (t1(i) <= 30)
    y3(i+1) = y3(i)/(1+h/2);
    t3(i+1) = t3(i) + h;
    i = i + 1;
end

%% Implicit solution for h=2
y4(1) = 1;
h = 4.2;
t4(1) = 0;
i = 1;
while (t4(i) <= 30)
    y4(i+1) = y4(i)/(1+h/2);
    t4(i+1) = t4(i) + h;
    i = i + 1;
end

%% Plot implicit solutions
figure(2); clf;
plot(t_true,y_true)
hold on
plot(t3,y3,'k.--','MarkerSize',20)
plot(t4,y4,'rs--','MarkerSize',10)
axis([0 20 -0.1 1.1])
xlabel('t')
ylabel('y(t)')
legend('Exact','Implicit Euler, h=1','Implicit Euler, h=4.2','Location','NorthEast')
