function[x,y]=odesolver(f,t0,t1,y0,steps)

h=(t1-t0)/steps;

x=[t0 zeros(1,steps)];
y=[y0 zeros(1,steps)];

for k=1:n
    x(k+1)=x(k)+h;
    y(k+1)=y(k)+h*f(x(k),y(k));
end

end