close all
clear all
T = 1; N=10; h = T/N;
t = linspace(0,1,N+1)';
f = @(t,y)-y.^2+2*t.*y;
y = @(t) exp(t.^2)./(1+sqrt(pi)*erfi(t)/2); %exact solution

yFE = zeros(size(t));yFE(1)=1; %forward Euler 
yBE = zeros(size(t));yBE(1)=1; %backward Euler

tFine = linspace(0,1,100)'; % fine mesh for plotting

for n=1:N
yFE(n+1) = yFE(n)+h*f(t(n), yFE(n)); 
yBE(n+1) =  (2*t(n+1)-1/h + sqrt((2*t(n+1)-1/h)^2+4*yBE(n)/h))/2; % see book for reasoning

plot(tFine,y(tFine),'b',t(1:n+1),yFE(1:n+1),'r-o', t(1:n+1), yBE(1:n+1),'k-x')
axis([0 1 0.75 1.18])
xlabel('t','fontsize',14)
ylabel('y','fontsize',14)
legend({'Exakt losning', 'Euler framat', 'Euler bakat'},'location','NorthWest','fontsize',14)
hold on
pause()
end


%% Add direction field 
[tMesh, yMesh] = meshgrid(0:0.1:1., 0.75:0.02:1.2);
F_t = ones(size(tMesh)); 
F_y = f(tMesh, yMesh);
figure(2)
quiver(tMesh, yMesh, F_t, F_y,1.,'b')
axis tight; xlabel('t'), ylabel('y')
title('Direction field for dy/dt = -y^2+2ty. That is (1,f(t,y))')
hold on

for n=1:N
plot(tFine,y(tFine),'b',t(1:n+1),yFE(1:n+1),'r-o', t(1:n+1), yBE(1:n+1),'k-x','linewidth',2)
axis([0 1 0.75 1.2])
xlabel('t','fontsize',14)
ylabel('y','fontsize',14)
legend({'(1,f(t,y))', 'Exakt losning', 'Euler framat', 'Euler bakat'},'location','NorthWest','fontsize',14)
hold on
pause()
end