function [x,y]=heun(f,a,b,ya,n)
h=(b-a)/n;
x=zeros(1,n+1);
y=zeros(1,n+1);
x=a:h:b;
y(1)=ya;
for i=1:n
    k1=f(x(i),y(i))
    k2=f(x(i)+h,y(i)+k1*h)
    y(i+1)=y(i)+h*(k1+k2)/2;
end