You learn from calculus that the derivative of a smooth function f(x), defined on some interval (a,b), is another function defined by the limit (if it exists) 

function H=heun(f,a,b,ya,m)
% Input -- f is the slope function entered as a string 'f'
%       -- a and b are the left and right end points   
%       -- ya is the initial condition y(a) 
%       -- m is the number of steps 
% Output-- H=[X',Y'] where X is the vector of abscissas and 
%          Y is the vector of ordinates 
h=(b-a)/m; 
X=zeros(1,m+1);
Y=zeros(1,m+1);
X=a:h:b;
Y(1)=ya;
for j=1:m
    k1=feval(f,X(j),Y(j));
    k2=feval(f,X(j+1),Y(j)+h*k1);
    Y(j+1)=Y(j)+h*(D(1)+(h/2)*(k1+k2);
end
H=[X',Y'];