function [x, y] = improvedEuler(fname, inte, y0, h)
x=inte(1):h:inte(2);
y(1)=y0;
for n=1:length(x)-1
y(n)=vpa(y(n),6);
k1=feval(fname,x(n),y(n));
y(n+1)=y(n)+h*k1;
k2=feval(fname,x(n+1),y(n+1));
y(n+1)=y(n)+h*(k1+k2)/2;
end
x=x'
y=y'
end
format long;
fanme=inline('2/3*x/(y*y)')
inte=[0,1];
y0=1;
h=0.1;
[x,y]=improvedEuler(fanme, inte, y0, h);
y=roundn(y,-5)
disp([x,y])
结果