% Euler's and Midpoint Method script
% solve IVP: 
% dy/dt = f(t,y), y(a) = y0

tic % timing your code, start of timing
%close all; clear all; % clear workspace, close open figures
N = 100; %number of steps
a = 1; b = 10; y0 = 0; h = (b-a)/N; %timestep
%f = @(t,y) (2./t).*y + (t.^2).*exp(t); %fctn handle to RHS of DE 
%ytrue = @(t) (t.^2).*(exp(t) - exp(1));
f = @(t,y) (-1/2)*exp(t/2).*sin(5*t) + 5*exp(t/2).*cos(5*t) + y;
ytrue = @(t) exp(t/2).*sin(5*t);

y = zeros(N+1,1);
t = [a:h:b]'; %turn into a column vector for fun!
y(1) = ytrue(a); %initial value, vector indexing. y is not a function handle

%Euler's method update
for i = 1:N
    y(i+1) = y(i) + h*f(t(i),y(i));
end

%midpoint method update
w = zeros(N+1,1);
t = [a:h:b]'; %turn into a column vector for fun!
w(1) = ytrue(a); %initial value, vector indexing. y is not a function handle

for i = 1:N
    w(i+1) = w(i) + h*f(t(i)+h/2, w(i)+(h/2)*f(t(i),w(i)));
end
toc % timing your code, end of timing

figure(1);
subplot(2,1,1);
plot(t,ytrue(t),'b-',t,y,'r-.'); title('Euler''s Method');
legend('true','approx','Location','northwest');
subplot(2,1,2);
plot(t,ytrue(t),'b-',t,w,'r-.'); title('Midpoint Method');
legend('true','approx','Location','northwest');






