f=@(x,y) -(y^2+x*y+x^2)/x^2;
% lahendame täpselt
% [x,y]=ode45(f,[x_0,x_n],y_0)
[x,yline]=ode45(f,[1,8],5) % saime, et y_n=-4.4365
ytapne=-4.4365;

n=1000;  % osalõikude arv
x0=1;
xn=8;
h=(xn-x0)/n; % sammude arv
% algtingimus
x(1)=1;
y(1)=5;

% Lahendame Euleri meetodiga
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)
x=[x0:h:xn];
plot(x,y,'o-')
% Lahendame prognoosi-korrelatsioonimeetodiga meetodiga
for i=2:n+1
  x(i)=x(i-1)+h;
  y(i)=y(i-1)+h/2*f(x(i-1),y(i-1))+h/2*f(x(i),y(i-1)+h*f(x(i-1),y(i-1)));
end
yprognoos=y(n+1)
hold on
plot(x,y,'*-')
% vordleme tulemusi [ytapne yeuler yprognoos]
[ytapne yeuler yprognoos] 