function hw1()
h1 = 0.19;
h2 = 0.21;

t1 = 1:h1:25;
t2 = 1:h2:25;

n = length(t1);
y1 = zeros(1,n);
y1(1) = 1;
for i = 2:n;
    y1(i) = y1(i-1)+h1*hw1func(t1(i-1),y1(i-1));
end

n = length(t2);
y2 = zeros(1,n);
y2(1) = 1;
for i = 2:n;
    y2(i) = y2(i-1)+h2*hw1func(t2(i-1),y2(i-1));
end

[T,Y] = ode45(@hw1func,[1 25],1);
plot(t1,y1,'r',t2,y2,'b',T,Y,'y');
legend('solution with h = 0.19','solution with h = 0.21','solution with ode solver');

function y = hw1func(t,x)
y = -5*t*x^2+5/t-1/(t^2);