%
%  Math 340 * Victor Matveev * April 20, 2012
%  See solution to homework #10 for details
%
%  dY
%  -- = x - Y
%  dx
%
clear;
set(gca,'fontsize',12);

y0 = 0.5;
T = 0.6;

n1 = 3;
n2 = n1*2;

h1 = T / n1;
h2 = T / n2;

f = @(x,y)  x-y;
y1(1) = y0;
y2(1) = y0;


for i = 1 : n1
    x = (i-1) * h1;
    y1(i+1) = y1(i) + h1*f(x,y1(i));
end

for i = 1 : n2
    x = (i-1) * h2;
    y2(i+1) = y2(i) + h2*f(x,y2(i));
end


x1 = h1 * [0 : n1];
plot(x1,y1,'ok--', 'markerfacecolor', [0 0 0], 'markersize', 6);
hold on;

x2 = h2 * [0 : n2];
plot(x2,y2,'db--', 'markerfacecolor', [0 0 1], 'markersize', 6);

yR = 2*y2([1:2:n2+1])-y1;
plot(x1,yR,'sk', 'markerfacecolor', [1 0 0], 'markersize', 6);

xx = linspace(0,T,100);
yy = (y0+1)*exp(-xx)+xx-1;
plot(xx, yy, 'm-', 'linewidth', 2);

legend('Euler h=0.2', 'Euler h=0.1', 'Richardson extrapol.', 'Exact');
xlabel('x'); ylabel('Y');
title('Solving dY/dx = x - Y');