dy = @(t,y)(-t^(-2) + 10*(y-1/t));
sol = @(t)(1./t);
h = 1;
t1 = 1:h:4;
y1 = zeros(length(t1),1);
y1(1) = 1;
for k = 2:length(t1)
    t = t1(k);
    ynew = ( y1(k-1) - h/t^2 - 10*h/t )/(1-10*h);
    y1(k) = y1(k-1) + h*dy(t,ynew);
end
plot(t1,y1,t1,sol(t1));
    
