% Ylesanded 16, ylesanne 1
f=@(x,y) -(y^2+x*y+x^2)/x^2;
[x,y]=ode45(f,[1,8],5)
plot(x,y,'o')
ytapne=interp1(x,y,8,'linear')

% a) Euleri meetodiga
x0=1;
xn=8;
y0=5;
h=0.1;
n=(xn-x0)/h;
x(1)=1;
y(1)=5;
y1=y0+h*f(x0,y0)

for i=2:(n+1)
   x(i)=x(i-1)+h;
   y(i)=y(i-1)+h*f(x(i-1),y(i-1));
end
x;
y;
yeuler=y(n+1)
hold on
plot(x,y,'--')

% b) Heuni meetod
x1=x0+h;
y1h=y0+h/2*(f(x0,y0)+f(x1,y1))

for i=2:(n+1)   
   y(i)=y(i-1)+h/2*(f(x(i-1),y(i-1))+f(x(i),y(i-1)+h*f(x(i-1),y(i-1))));
end
x;
y;
yheun=y(n+1)

hold on
plot(x,y,'*')

% c) keskpunkti meetod
y1k=y0+h*f(x0+h/2,y0+h/2*f(x0,y0))

for i=2:(n+1) 
   y(i)=y(i-1)+h*f(x(i-1)+h/2,y(i-1)+h/2*f(x(i-1),y(i-1)));
end
x;
y;
ykeskp=y(n+1)

hold on
plot(x,y,'diamond')

[ytapne yeuler yheun ykeskp]